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I.  INTRODUCTION 


Numerical  weather  prediction  with  the  aid  of  primitive 
equations  has  now  been  performed  operationally  for  several 
years  (Reiser,  1969;  Shuman,  1968).  Forecasts  for  one  and 
two  days  by  primitive  equations  indicate  some  improvement 
when  compared  with  forecasts  with  the  quasi-geostrophic  and 
filtered  equations.  However,  it  has  not  been  shown  in  a 
clear  and  convincing  way  whether  this  is  solely  due  to  the 
unfiltered  part  of  the  equations  or  to  purely  numerical 
improvements  such  as  improved  vertical  resolution,  an 
alternating  grid,  and  the  introduction  of  new  physical  effects 
such  as  sensible  heat,  latent  heat,  radiation,  etc.  When  the 
forecasts  with  the  primitive  models  are  extended  further  in 
time,  the  improvements  achieved  with  these  models  seem  to  be 
more  obvious . 

If  a  quasi-geostrophic  model  is  compared  with  a  primitive 
model  with  a  resolution  of  three  to  four  vertical  layers  or 
less,  it  will  be  found  that  the  computational  time  (and  cost) 
for  the  forecast  with  the  primitive  model  is  roughly  about 
ten  times  as  large  as  for  the  quasi-geostrophic  model,  if 
explicit  time-integration  schemes  are  used.  The  grid  lengths 
now  in  use  for  operational  weather  prediction  are  still  about 
300  to  400  km. 


3 


Since  significant  weather  disturbances  have  dimensions 
which  are  only  three  to  five  times  that  size,  it  is  obvious 
that  the  truncation  errors  in  the  computation  of  the  horizon¬ 
tal  finite  differences  are  much  too  large. 

That  grid  distances  of  300  km  and  more  have  been  used 
for  such  a  long  time  is  naturally  due  to  insufficient  com¬ 
putational  capacity,  but  may  also  partly  be  due  to  accustomed 
routine.  However,  a  decrease  in  the  horizontal  grid  length 
to  half  the  size  implies  an  increase  in  the  number  of  grid 
points  by  a  factor  of  four  for  the  same  computation  area. 

Due  to  the  criterion  of  Courant,  Friedrichs  and  Lewy , 
the  time  step  must  be  decreased  by  a  factor  of  two  in  order 
to  maintain  computational  stability.  If  the  computations 
can  be  organized  in  the  same  way  as  for  the  larger  grid,  a 
halving  of  the  grid  length,  therefore,  implies  an  eight-fold 
increase  in  the  computation  time.  From  the  operational  point 
of  view,  therefore,  a  primitive  model  should  be  compared  with 
a  quasi-geostrophic  model  where  the  horizontal  grid  length 
has  been  decreased  to  half  the  size. 

When  the  forecasting  areas  are  small  and  cover  only  one 
third  or  less  of  a  hemisphere,  the  horizontal  boundaries  will 
fall  in  meteorologically  active  areas.  This  disadvantage  does 
not  create  any  large  problem  for  the  filtered  models  and  a 
moderate  horizontal  smoothing  in  the  neighborhood  of  the 
boundaries  is  sufficient.  For  the  primitive  models  the 
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problem  is  much  worse ,  since  high-amplitude  gravity  waves  are 
generated  at  the  boundaries.  These  waves  propagate  into  the 
area  and  greatly  affect  the  meteorological  information. 

Except  for  some  successful  experiments  (Bushby,  1967,  1968; 
Gerrity,  19  69)  ,  there  has  not  been  reported  any  adequate 
technique  to  avoid  this  in  the  general  case.  For  this  reason, 
forecasts  for  restricted  areas  with  the  complete  equations  in 
operational  use  imply  considerable  difficulties. 

It  may  now  be  argued  that  it  is  of  no  use  to  apply  a 
quasi-geos trophic  model  to  a  fine  mesh,  since  that  means  the 
model  will  predict  (or  try  to  predict)  scales  of  motion 
characterized  by  large  Rossby  numbers.  For  instance,  when 
the  Rossby  number  is  on  the  order  of  one,  all  the  terms  in 
the  vorticity  equation  are  of  equal  magnitude.  The  quasi- 
geostrophic  models  will  thus,  according  to  this  analysis, 
give  rise  to  intolerably  large  errors  for  small  and  intense 
vortices,  especially  at  low  latitudes.  However,  it  has  been 
shown  (Bengtsson  and  Moen,  1971)  that  substantial  improvements 
in  forecasts  from  quasi-geostrophic  models  are  obtained  if 
the  grid  size  is  reduced  from  300  to  150  km. 

The  reason  for  this  is  that  the  higher  order  terms  in 
the  vorticity  equation  to  a  considerable  degree  cancel  each 
other.  It  is  only  in  the  final  stage  of  the  cyclone  develop¬ 
ment,  when  the  flow  becomes  very  deformed,  that  these  terms 
become  important. 
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It  is  also  possible,  as  will  be  shown  in  this  report,  to 
include  these  terms  in  the  integration  and  estimate  them  with 
the  aid  of  quasi-geos trophic  divergence. 

It  is  the  experience  of  the  author  that  filtered  models 
still  are  very  useful  for  short-range  predictions  at  medium 
and  high  latitudes.  It  is  also  very  probable  that  unsuccess¬ 
ful  predictions  of  especially  rapid  cyclogenesis  are  due  to 
inaccuracies  in  the  initial  state  and  to  unsatisfactory  ways 
of  including  topographical  effects ,  parameterization  of 
dissipation,  and  the  heating  mechanisms.  Recent  comparisons 
performed  in  Sweden  between  filtered  and  primitive  equation 
models  support  this  view. 
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2.  PROGNOSTIC  EQUATIONS 


The  vorticity  equation,  thermodynamical  equation  and 
the  continuity  equation  read: 


3£  = 
3t 


-  WVrj  -  f  D; 


a  ,  a  x  iTrv7/^^\  n  /  p  R  ,<sqx  . 

It  (I?>  =  -''V-7(I?)  ■  R(rd"r)p  ■  S3  (dt); 


3w 

Ip 


4rr  =  -  D; 


where 


\V  =  lkXV<j> , 

D  =  V*\V, 

r ,  =  — —  ,  and 
d  pc 


r  _  8T 

ap 


(2.1a) 

(2.1b) 

(2.1c) 


We  will  now  introduce  some  special  model  assumptions.  We 
assume  that  the  atmosphere  is  bounded  by  a  pressure  surface 
near  the  surface  of  the  earth,  p=pQ,  and  an  upper  pressure 
surface,  p=p^,  near  the  tropopause.  We  will  further  assume  a 
third  interjacent  pressure  surface,  p=pm,  which  separates  the 
atmosphere  in  two  layers.  We  will  now  represent  the  wind 
field  in  the  following  way: 
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\V  =  XV  -  2\VX 

^  o  x  m 

V  =  'vm  +  »2 

Mn  rl 


\V  =(Wm  +  2V2)  2_ 


layer  1 
layer  2 
layer  3 


Figure  1.  Vertical  representation  of  the  wind  for  a  3-para- 
meter  model  . 

According  to  the  definition,  equations  for  D  and  £  will  have 
the  same  form  as  those  for  V.  Equation  (2.1a)  can  now  be 
integrated  between  P=PQ  und  p=0  with  the  boundary  conditions 
w(p=0)=0  and  w(p  )=w  . 

This  gives  a  prognostic  equation  for  the  vertically 
integrated  vorticity: 


It  (Cm  +  °1?2  -  c251>  =  -  VV(c3rVc2?l+c4C2+c7f) 

-  \V1-V(-c2nm+c5£1)  -  W2*V(c4nm+c6?2+2c7f)  +  Cgf03s  (2.2) 
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where 


2p. 


m 


°i  = 


6pm'2pl 
'4  SPo-JPj 


C7  6P0-4Pl 


= 


2 (p  -p  ) 

*0 

2pcTpl 


„  _  8po'8pm 

5  6po'3pl 


= 


•S  2po-Pl 


_  6pp-4pl 

3  6po"3pl 


8p 


= 


m 


6  6pQ-3p1 


The  next  two  prognostic  equations  are  computed  from  the 

difference  between  the  vorticity  equation  at  level  pm  and  pQ 

and  the  corresponding  difference  for  p,  and  p  .  We  will  now 

1  rm 

get  two  vorticity  equations  valid  for  layer  1  and  layer  2: 


0 V-\V  )-V?  +\V  •V(n  -£,)  +  f 

ml  l  ±  ml 


D, 


0; 


(2.3) 


Vt2  +  wm+v2)-vc2  +\V2 


v  (n  +C,)  +  f  D0  =  0 . 

m2  2 


(2.4) 


From  these  two  equations  we  will  now  eliminate  the  diver¬ 
gencies  and  with  the  aid  of  the  continuity  equation 
(2.1c)  and  the  thermodynamical  equation  (2.1b).  We  will  first 
integrate  the  continuity  equation  between  the  two  levels  p 
and  p^ : 


^b 

w  (p  )  =  w(pb)  +  f  D  dp. 

Pa 


(2.5) 
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We  will  now  put  p  =p ,  p,=p  and  w(p  )  =w  and  will,  thereby, 

3,  D  O  U  b 

get  an  expression  for  co  (p)  in  layer  1: 


co  (p) 


CO  + 

s 


(p0-p) 


(p-p  )  -  (p  -p  )  2 

^  ,  m  *0  cm  ^ 

D  +  -  D,  . 

TT1  p  — p  I 

ro  ^m 


(2.6) 


With  p  =p,  p  =p  and  co(p  )  from  (2.6)  we  get  the  following 
^a  b  ^m  m 

expression  for  co  in  layer  2: 


(Pm-P) 2 

“(P)  =  ws  +  (po_p)  Dm  "  (po"Pm)  D1  +  P~-P~  °2* 
With  pa=p,  Pb=0  and  co  (p=0 )  =0  we  will  have  for  layer  3: 

“<p)  *  '2^;  (Dm  +  2D2>  ' 


(2.7) 


(2.8) 


Finally  we  get  a  relation  between  D^,  ,  D ^  and  cog  with  the 

aid  of  an  integration  of  the  continuity  equation  from  the 
top  to  the  bottom  of  the  atmosphere: 


D  = 
m 


C2D1 


°1D2  - 


CoC0 

8  s 


(2.9) 


We  now  introduce  the  stream  function  into  the  thermo¬ 
dynamical  equation  and  integrate  through  layers  1  and  2 . 

(T  -r)  is  assumed  to  be  constant  in  every  layer  and  (^) 

1:1  pm 


6Q 

dt 


1 


0  . 
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dip. 


2f 


3t 


•2E'Vvh  +  R(r4-D1  ;°^  dp  +  ln&  <|£)  ; 


m 


p 
p 


2f 


9^: 

3t 


m 


■2£Vm-V<<2  +  E(rd-D2  / 


CO 


dp 


(2.10) 


The  integrals  /  —  dp  can  be  computed  with  the  aid  of 

P 

equations  (2.6),  2.7)  and  (2.9),  and  the  system  (2.10)  can  be 
written 

mlDl  +  m2°2  =  Hl; 


nl°l  +  n2°2  =  H2; 


(2.11) 


where 


mi  =  I  (rd-r)i 


2  2 

p.  (p  -p  )  +p  (2p  -p. )  p 

1  *o  *m  ^m  v  ^o  ^1'  ,  ,o . 

-  •  In  ( — ) 
P 


(p  -p  ) (2p  -p. ) 

^o  Mn  ^o  *± 


2 ‘Pp-Pm1  1  ,  .  ,  , 

2po-pd  2  po  pm  pm 


itu  = 


T  <rd-r>i 


2p  p 

^o^m 

2p0‘pi 


p  2p  (p  -p  ) 

In  (— )  +  mom 


m 


2po-pl 


ru  =  ■ 


f  «d-r)2 


p.  (p  -p  )  p  2  (p  -p  )  (p  -p  ) 

1  o  ^m  ,  ,  m.  o  ^m  ^m 

In  (—  -  — 


2p  -p. 

o  1 


2po~pl 


1 1 


ru  - 


?  (rd-r>2 


ElV.fg°3;!  ln  ,V 

(2p0-Pl)  (Pm-Pl>  Pl'  2po-Pl 


+  (p  +p,  )  -  2p 

2  ^m  rl 


dip 

H1  =  £  ST 


pm  dt  P 
m  •Lo 


i  <rd-r)i 


2Po'Pl 


p  2 (p  -p  ) 

ln(-^)  +  °  m 


m 


2Po‘Pl 


COc 


dip. 


R 


H2  =  f  9t  +  fVV*2  -f  (rd-r,2 


P 

In  (-Si) 


2po-Pl  Pl 


+ 


2(pm-pl> 

2£vpt 


CO, 


From  the  system  (2. 11) we  can  easily  express  the  diver¬ 
gencies  in  an  explicit  way: 


D1  " 


= 


“alHl  +  a2H2; 


blHl  b2H2; 


(2.12) 


where 


-n. 


-nr 


a..  = 


b.  = 


mln2"nlm2 


-n. 


mln2"nlm2 


a_  - 


b„  = 


mln2"nlm2 


-m. 


mln2  nlm2 


12 


Introducing  these  expressions  for  and  D2  into  the 
prognostic  equations  (2.3)  and  (2.4)  and  expressing  the  wind 
and  vorticity  in  terms  of  the  stream  function  gives: 


„2  9*1  a  f2  +  a  f2  f*2 

v  at  if  'at  2f  at 


-  -  “  JWi;nm”2?l) 


+  a1f2J(^m;ip1)  -  a2f2j^m;4'2)  "  a3fu)S  ~  alfH  (2.13) 


2  ^2  2  3^2  2  3*1 
»2  UT  -  b2f  TPT  +  blf  — 


=  -  J<VC2>  •  J<','2;nm+2i:2, 


-  b1f2J(^m;^1)  +  b2f  J(^m;^2)  +  b3fu)s  +  blfH  (2.14) 

Here  we  have 


a3  alSl  a2s2; 


b3  blSl  b2S2; 


where 


R, 


Si  =  flVD! 


S2  =  f (rd"r) 2 


,  p  2(p  -p  ) 

1  In  (— )  +  °  m 


2pcTpl  pm 
P 


2pcTpl 


Pm  2 (p  -p,) 
ln(— )  + 


2po~pl  pl  2po-pl 


H  =  in  (_°)  (iS)  =  h.  (|2)„  with  h  -  ln(— )  . 

4cp  pm  dt  po  1  dt  po  1  4cp  pm 

We  will  also  introduce  the  stream  functions  into  equation 

(2.2)  and  define  by  ^M=^m  ~  c2^l+c1^2’ 
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This  results  in  the  equation 
2 

’5T  =  -  J(Vc3nnfc2?l+c4i;2+c7£)  '  J<*1'_c2,Vf<V-l) 

-  J(iK2;c4nm+c6;2+2c7f)  +  c^.  (2.15) 

The  equations  (2.13)  through  (2.15)  now  constitute  our 

system  of  prognostic  equations.  The  only  thing  which  we  now 

have  to  do  is  to  find  an  expression  for  the  non-adiabatic 

heat  H,  and  the  lower  boundary  condition,  w  . 

s 
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3.  BOUNDARY  CONDITIONS 


3.1  LOWER  BOUNDARY  CONDITIONS 

We  now  assume  that  can  be  separated  into  two  parts; 

one  part  which  depends  upon  the  dissipation  in  the  boundary 

layer  w'  and  one  part  which  depends  upon  topography  a)".  We 

s 

thereby  assume: 

us  =  “s  +  w"  .  (3.1.1) 

From  the  Ekman  theory  about  the  variation  of  the  wind  in 
the  friction  layer,  we  could  easily  derive  the  following 
expression : 

=  -  gP0  y|T-F.C0  where  F  =  1  +  c*sin0  -  c*cos0, 

the  wind  in  the  surface  layer  has  a  magnitude  c |w  | /  the 
angle  between  the  geostrophic  wind  and  the  surface  wind  is 
given  by  0,  K  is  the  turbulent  coefficient  of  the  viscosity, 
and  po  is  the  densitY  of  the  air  at  the  surface.  Inserting 
the  following  numerical  values : 

K  =  10  m/s;  g  =  9.81  m/s2;  f  =  10_4s-1;  Tq  =  280°K; 

R  =  287;  pQ  =  100  cb;  gives 

=  “  2.729597  F*Cq  =  -2.729597^(^-2^).  (3.1.2) 
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F(c,0)  can  be  given  a  constant  value  in  the  model  or  we  can 

also  assume  different  values  over  land  and  sea. 

The  following  values  will  be  used: 

Over  land  0  =  10°  and  c  =  0.78  gives  F  =  0.36635, 

that  is  o)'  =  -1.?  . 

s  o 

Over  sea  0  =5°;  c  =  0.85  gives  F  =  0.22734, 


that  is.o)'  =  0.62055  £ 

'  s  o 

If  is  assumed  to  be  that  part  of  the  air  which  is 


forced  over  the  mountains ,  we  get 


oo 


k2  V’Ps  =  k2  'VP£ 


(3.1.3) 


=  1  will  be  used  in  the  model. 


The  equations  (3.1.1  -  3.1.3)  now  give  the  lower  boundary 
condition  for  w  .  For  further  information  see  Bengtsson 
(1969) .  This  way  of  treating  the  topographical  effect  as  given 
by  equation  3.1.3  is  quite  unrealistic  and  seems  to  underesti¬ 
mate  the  effect  of  the  topography.  (This  will  be  especially 
true  for  steep  and/or  small  scale  mountains.)  A  new  way  to 
treat  mountains  as  impenetrable  vertical  barriers  has  recently 
been  published  (Egger,  1972) .  A  similar  way  to  include 
mountains  in  vertically  integrated  balanced  models  will  be 
described  by  the  author  in  a  coming  investigation. 


3.2  LATERAL  BOUNDARY  CONDITIONS 

The  model  can  use  two  different  kinds  of  lateral  boundary 
conditions  namely : 
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a. 


Constant  Inflow 


_11  =  _2  =  _JH  =  o 

3 1  3 1  3 1 

^(t)  =  C1(t=0)  (=  constant  in  time) 

£2(t)  =  ?2(t=0)  (=  constant  in  time)  (3.2.1) 

£  (t)  =  t,  (t=0)  (=  constant  in  time) 

m  m 

b.  Variable  Inflow 

3^m 

The  values  for  ,  5^,  £2 ,  and  along 

the  boundary  are  generated  through  an  integration  for  a 
larger  area  which  includes  the  actual  area.  Interpolation 
in  time  and  space  is  necessary  to  syncronize  the  boundary 
values.  A  technical  description  of  this  is  found  in  section 
10. 
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4. 


PARAMETERIZATION  OF  PHYSICAL  PROCESSES 


4.1  SENSIBLE  HEAT 

In  the  computation  of  equation (2 . 10)  we  assume  that  ^ 
decreased  linearly  to  0  from  pQ  to  pm.  The  sensible  heat 
which  is  transported  to  the  atmosphere  from  the  underlying 
surface,  is  introduced  in  the  following  way: 

(4?)  =  ( Aq | v0 | +A2 )  (Ts-To)  over  sea,  if  Ts>Tq; 

po 

(4.1.1) 

(^Q)  =  o  over  land  and  over  ocean  areas  if  T  <T  ; 

dt  p  so 

^o 

where  Aq  and  A 2  are  empirical  constants  A2=10Aq 

_  9  -l 

=  5.10  m (sec)  (deg)  .  Tg  is  the  sea  surface  tempera¬ 
ture  and  Tq  is  the  air  temperature  near  the  sea  surface. 

An  approximative  temperature  at  the  level  Px  =  j  (pQ+PM) 
is  obtained  through  an  integration  of  the  hydrostatic  equation: 


2f 


T  = 
x 


R  in  (^) 

^m 


(4.1.2) 


We  now  assume  that  T  is  constant  in  the  layer  and  is  75%  of 
rd.  We  therefore  get 


2f 


To  =  Tx+0'75  rd  = 


po 

R  In  (— ) 
P 


1  + 


°.?5  R(Pc-Pm) 


c  (p  +P  ) 

p  ^o  Mil 


1 . 


(4.1.3) 


I  8 


In  equations  (2.13)  and  (2. 14) the  sensible  heating  (H)  was 
defined  to  be: 


H  = 


R 


4C 


,  ,Po,  ,SQ> 
ln  <p5>  <at> 


Po 


Substitution  of 
expression  for  T 


(^r)  from  equation  (4.1.1)  and  using  the 
Po 

from  equation  (4.1.3)  gives: 


H  =  h  -10  2  (0 . 5  •  10  1|V  1-0.5)  (T  -h  '\p  ) 

X  O  S  Z  1 


where 


,  R  ,  ,  o, 

hl  =  4V  ln  (f’  ! 

p  ^m 


2f 


h„  = 


R  in  (^) 
^m 


1  + 


0.75  R (p  -p  ) 
_ ro  ^m 

c  (p  +p  ) 
p  o  ^m 


(4.1.4) 


4.2  PROGNOSTIC  EQUATION  FOR  HUMIDITY 

The  specific  humidity  can  be  predicted  by  the  following 
equation: 

-|SL  =  -  \Wq  -  +  e  ~  r  +  AV2q  .  (4.2.1) 

Here  e  denotes  evaporation,  r  condensation,  and  A  is  the 
coefficient  of  dissipation.  We  will  disregard  e  and  r  will 
be  introduced  in  a  different  way.  Using  the  continuity 
equation  we  get  the  following  prognostic  equation: 

-|2-  =  -  V*  (q  \V)  -  |p  (qod)  +  AV2q.  (4.2.2) 
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Since  our  model  has  a  very  low  vertical  resolution  we  have  to 
parameterize  the  vertical  distribution  of  humidity.  We  will, 
therefore,  use  precipitable  water  as  the  prognostic  variable. 
The  precipitable  water  is  defined  as 

p 

00  ^o 

p  =  /  p  dz  =  —  /  q  dp  (4.2.3) 

^w  Fw  g  ^  r 


where 

pw  =  the  density  of  water  vapor  and  pT  is  the  pressure  at 
the  level  over  which  we  can  disregard  the  humidity.  (Here 
we  have  p^  =  30  cb  and  pQ  =  100  cb . ) 

A  new  quantity 


w  = 


P0'PT 


(4.2.4) 


which  may  be  called  normalized  precipitable  water  is  intro¬ 
duced.  We  assume  that  q(x,y,p,t)  =  gE (p)w(x,y ,t)  where  E(p) 
describes  the  vertical  variation  of  q  computed  from  the 
standard  atmosphere  under  the  assumption  that  the  relative 
humidity  is  50%.  The  expression  (4.2.4)  is  now  introduced 
into  (4.2.2)  and  we  then  integrate  with  respect  to  p  from  pQ 
to  to  give 


^  =  -  V*(\Vw)  -  wd'w  +  AV2w 

d  t  S 

where 

1  p° 

y  =  - -  /  $7(p)  E(p)  dp  and  d' 

Po'PT  pT 


E<Pq> 

po_PT  ' 


(4.2.5) 


20 


With  \V  =  IkxVip  +  Vx  equation  (4.2.5)  can  be  written 

2  2 

—  =  -  J($,w)  -  V x •  Vw  -  wV  X  -  wd'ii)  +  AVw 

9 1  s 


(4.2.6) 


where 


ip  = 


e  \p  +  e-jip-,  + 
nrm  11 


e2*2 


and 


V  X 


=  D  = 


e  D  +e1D1+e0D„. 
mm  11  2  2 


4 . 3  LATENT  HEAT 

There  are  two  conditions  for  condensation: 

(a)  <  wtQl  (where  w  is  a  given  tolerance) 


(b)  the  relative  humidity  should  exceed  and  be  equal 
to  80%. 

If  these  two  conditions  are  valid,  the  latent  heat  is 
computed  by  the  aid  of  expression  (4.3.1) 

P, 


H 


-  —  In  (— )  (— )  -  h  (— ) 

lat  4cp  xn  dt  lat  nl  Wlat 


where 


(3t>lat  =  '  L'“lF  and  F  = 


X 


F  +  §-«l? ) 

p  p 


~  X 

t  (ia_) 

U3p  JT 


(4.3.1) 


+ 


R  T  ,9q*  , 

Cp  p  '  3T  V 


X 

q  is  the  maximum  specific  humidity. 

and  (b)  are  not  valid,  we  put  H.  = 

X  St 


If  the  conditions  (a) 
0.  Also  see  paragraph 


8.4. 
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5.  COMPUTATION  OF  THE  VERTICAL  MOTION 


The  integrated  vertical  motions  in  the  two  layers  1  and 
2  are  computed  in  the  model.  With  the  aid  of  the  equation 
(2.6)  and  (2.9)  we  obtain: 


03.  = 


1  to 


Pn  P  +Pm"Pn  ■=>•(  2P-,  -3p  -p  )  (p  -p  ) 

f°  codp  =  . °  m  1  03  +  3 _ -1 _  m  o  p 


1  Po"Pm  pm  '  2Po‘Pl “  ”s  ■  2po-P! 

Pm(Po-Pm>  „ 

2p0-P1  D2 1 


(5.1) 


“l  =  tl“s  +  t2Dl  -  *3*2' 


In  the  same  way  we  obtain  from  the  equation  (2.7)  and  (2.9): 


W  =  - 

2  P  “P 


P  P 

tm  ,  ^m 

J  03dp  = 


03 


m  rl  p. 


2p  -p  s 
o  1 


P  (p  -p  ) 
m  *0  ^m 

2p  -p. 

ro  *1 


D1  + 


+ 


3(2Po-Pl>  ‘Pm-Pl’-Pm^Po-Pm-Pl1 


2p  “Pi 

o  1 


V 


(5.2) 


032  =  t403s  -  t3D1  +  t5D2  . 


Here  we  have 

P  +P  -p, 

.  _  o  ^m  ^1 

1  2p  -p 

o  rl 


t  =  i  (2pr3pm~Po)  (p0-pm> 

23  2p  -p 

o  ^1 


t  _  Pm^o-Pm1 

3  2po‘pl 
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t  =  -m— 

4  2P0-Pl 

.  _  |(2pQ-Pl)  (pm-Pl)  -  Prn(2po-pm-p1) 

5  (2P0"P1) 


The  physical  parameters  are  computed  by  the  subroutine  COEFF3P. 
(See  appendix  B  to  this  report) 
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6.  NUMERICAL  VALUES  OF  THE  CONSTANTS 


For  the  levels  pQ  =  100  cb ,  Pm  =  50  cb ,  p  =  30  cb  and 
the  stabilities  (T^-r^  =  0.422222,  (rd~r)2  =  0.511111  we  get 

the  following  numerical  values  for  the  constants: 


c1  =  yj  =  0.588235  ir^ 

c2  =  —  =  0.588235  m2 

c3  =  =  0.941176  n1 

c4  =  =  0.470588  n2 

40 

c.  =  ~  =  0.784314  a. 

5  51  1 

cr  =  ~  =  0.784314  a„ 

6  61  2 


-  826.330 

-  688.40 

-  532.92 

-  10.58.36 

2.0828*10" 

1.3546*10" 


c7  =  -gy  =  0.588235  b, 
C8  =  I70  =  °*0117647  b2 


1.0475*10 

1.6261*10 


s1  =  28.230856 
s  2  =  10.645832 

a3  =  0.044374 
b3  =  0.012258 

h±  =  0.495351*10 

h2  =  0.110953*10' 

h3  =  1.03552*10"^ 
h4  =  0 

h5  =  0.492929*10“ 

hr  =  1.03552*10"^ 
t,  =  0.705882 

t2  =  -18.627500 

t3  =  14.705900 

t4  =  0.294118 

tc  =  -28.627500 
5 
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7. 


INTEGRATION  OF  THE  COMPLETE  VORTICITY  E&UAT ION 


7 . 1  GENERAL  ASPECTS 

Very  little  knowledge  exists  about  the  effect  of  the 
small  order  terms  in  the  vorticity  equation:  the  advection 
of  vorticity  by  the  divergent  wind,  \Vx*V£;  the  product  of 
relative  vorticity  and  divergence,  £V*  \V;  the  vertical 
advection  of  vorticity,  and  the  twisting  term  Ik •  (-!— jxVw)  . 

If  these  terms  are  used  it  is  necessary,  in  order  to 
conserve  the  total  energy  for  an  adiabatic  model,  to  use  the 
complete  balance  equation.  If  this  is  not  the  case,  the  model 
will  not  conserve  total  energy  and  after  a  certain  time  the 
development  starts  to  deteriorate.  This  judgment  has  been 
mainly  qualitative  and  we  do  not  know  the  size  of  the  error 
due  to  this  inconsistency.  It  may  be  that  this  error  is 
relatively  small  in  comparison  to  other  errors,  as  for  instance, 
uncertainties  of  the  initial  state,  and  uncertainties  in  the 
description  of  the  dissipation  and  the  heating  mechanisms. 
Therefore,  it  is  necessary  to  perform  a  more  detailed  study 
of  the  problem  and  base  our  decision  on  a  quantitative 
investigation . 

It  is  by  no  means  evident  that  we  should  use  the  same 
kind  of  assumptions  in  the  formulation  of  models  for  short- 
range  predictions  (24  hours)  as  for  models  for  medium-  and 
long-range  prediction.  An  example  will  illustrate  this. 
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If  one  is  interested  in  long-time  integrations  of  the  baro- 
tropic  vorticity  equation,  it  is  necessary  to  use  a  finite 
difference  expression  which  conserves  kinetic  energy, 
vorticity,  and  mean-squared  vorticity.  A  finite  difference 
expression  which  conserves  these  identities  has  been  derived 
by  Arakawa.  However,  if  one  is  interested  in  short-range 
predictions,  the  so-called  "Arakawa  Jacobian"  is  not 
recommended  since  the  phase-speed  error  is  larger  than  the 
conventional  finite  difference  analog  to  the  Jacobian  operator 
which  only  conserves  vorticity.  Experiments  have  shown  that 
for  forecasts  up  to  four  or  five  days  it  is  not  necessary  to 
use  an  energy  consistent  Jacobian  operator  since  the  error  in 
the  kinetic  energy  is  much  smaller  than  errors  due  to  other 
effects . 

One  of  the  problems  which  we  have  with  the  simplified 
vorticity  equation  is  the  over-prediction  of  anticyclogenesis. 
This  seems  due  mainly  to  the  lack  of  the  term  £ V *V  in  the 
vorticity  equation: 

|§  =  -  -(f+c) V-W.  (7.1.1) 

In  areas  of  convergence,  relative  vorticity  is  mostly  positive, 
or  will  be  after  a  short  time.  This  means  that  the  relative 
vorticity  will  increase  faster  in  such  areas  if  the  complete 
expression  (7.1.1)  is  included.  On  the  other  hand,  in  areas 
of  divergence,  the  relative  vorticity  is  mostly  negative,  or 
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will  be  after  some  hours.  If  the  complete  expression  is  used, 
the  relative  vorticity  will  decrease  more  slowly  than  if  we 
use  the  simplified  expression.  It  is  easily  seen  that  this 
term  will  create  an  asymmetry  in  the  vorticity  pattern  which 
is  also  observed  in  reality. 

Also  the  vertical  advection  of  vorticity  seems  to  play 
an  important  role,  especially  in  cyclone  development.  During 
the  development  of  the  cyclone  the  activity  is  mainly  con¬ 
centrated  in  two  different  areas. 

One  area  is  in  front  of  the  warm  front  or,  later  in  the 
development,  the  occluded  part  of  the  front.  The  other  area 
is  found  below  the  upper-air  low  or  trough,  where  a  special 
center  of  activity  is  created  in  the  later  stages  of  the 
cyclone  development.  During  the  development  of  the  cyclone, 
an  area  of  sinking  motion  is  concentrated  under  the  upper- 
air  low. 

8?  _  _  3£ 

3t  W3p  (7.1.2) 

It  is  easily  seen  from  the  vorticity  equation  (7.1.2)  that 
this  effect  will  give  an  increase  in  the  relative  vorticity 
in  areas  of  sinking  motion  and  where  the  vorticity  increases 
with  height. 

Upward  motion  over  a  surface  low  yields,  in  the  same 
way,  an  increase  in  the  vorticity  for  the  levels  above. 
Therefore  the  vertical  advection  of  vorticity  will  increase 
the  speed  of  occlusion.  Figure  2  shows  two  different  12-hour 
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Figure  2.  Two  example  12-hour  predictions  with  a  5-layer  quasi-geostrophic  model. 


predictions  with  a  5-layer  guasi-geostrophic  model.  Solid 
lines  indicate  a  prediction  performed  with  a  quasi-geostrophic 

energy  consistent  model.  Dashed  lines  indicate  a  prediction 

3  £ 

where  the  terms  -(?V»  V  +  have  been  included. 


7.2  DERIVATION  OF  THE  FORECASTING  EQUATIONS 
The  complete  vorticity  equation  reads: 


3  V 


It  =  ~  "  \Vx’Vn  ■  fV*Wx  -  S  V*vx  -  0)||  +  Ik-(g^xVm).  (7.2.1) 


Here  V,  is  the  non-divergent  wind  and  V  the  divergent 

T  X  ' 

wind.  The  terms  1,  2,  3  and  4  are  denoted  non-geostrophic 
terms.  Integrating  through  layer  1  by  the  representation  of 
V,  Z,  and  D  given  in  section  2  yields 


(P~-P J  I 


3?  3?n 

1 


o  ^m  3 1 


3t 


=  'Po'Pm* 


+  wl‘v(sm+f> 


-  fwpV^  -  f(Dm-Dl)  -  ?m(Dm-Dl)  +  Cj.  (Dm  -  fop 

+2C1w1  -  2tk*(V1xVw1) 


Integration  through  layer  2  and  layer  3  yields  in  a 
similar  way 


(7.2.2) 


3?  3<So 

/  N  r  m  J  2 

Pm-Pl  3T- 


=  ‘Pm-Pl1 


-«V7(Wf>  -W2-v(Vf) 


|v2-v;2  -  f(Dm+D2)  -  ;m(Dm+D2)  -  c2<v!d2) 


+  2C2W2  “  2  ^2xVa)2^ 


(7.2.3) 
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2P1 


m  +  2-  2 


at 


at 


“  2P1 


-  (  Vm*AV2)-V£  -  j(  \Vm+2V2).V(Cm+2C2) 


-f(V2D2>  -  3<V2?2>  <V2D2> 


+  Ik*  i  Wm+2  \V2)  xVoo3 


(Cm+2?2)  W3 


(7.2.4) 


where  and  w2  are  given  in  (5.1)  and  (5.2)  and  is  given 
by : 


-  1  P1 
o)_  =  —  /  wdp 

3  p, 

^1  o 


(7.2.5) 


If  we  now  add  (7.2.2),  (7.2.3)  and  (7.2.4)  and  divide  by 

2p  -pn 

o  I  we  get  a  prognostic  equation  for  the  vertically  mte- 
2 

grated  mean  vorticity.  From  now  on  we  will  only  keep  the 
non-geostrophic  terms  on  the  right  hand  side  of  the  equation. 


at  (?m+cl?2~c2Cl)  *  (  VVX)m’V(c3T1nrc2Cl+c4?2+C7f) 

-(  l' v?(-c2Cm+C5’l)  "  (  VX*  2‘V(c4nm+C6?2  +  ^c7f)  *1 


+  {c8?ma3s+c2Cl(Dm  3D1}  C2  (c4Dm+C6D2)  +  ‘’m(c7Dm+2c7D2)  }2 


+c8{2C1w1+2S2w2  -  (<;m+2?2)  w3>3  +Cg{  (k*  [-2^1xVw1-2V2xVw2 

+  (NVm+2  \V2)xS]}4  (7.2.6 

The  two  remaining  prognostic  equations  are  computed  from  the 
difference  between  the  vorticity  equation  for  pm  and  pQ  and 
the  difference  between  the  vorticity  equation  for  p^  and  p^  . 
We  will  then  have  two  vorticity  equations  valid  for  layer  1 
and  layer  2  respectively.  The  geos trophic  terms  are  omitted. 
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3?- 


3^  =  -  {^X)m-(^X)l),VCl  +  (VX}rV(T1m-?l)}l 


{(Cm  2?1)  D1+C1Dm}2  +  {p  -p  (%  Us)}3 

o  m 


(7.2.7) 


3?. 


-  <<Wx)m+Wx>2>-V?2  + 


{  <Cm+2C2)  D2  +  i’2Dm)2  +  (p  -p,i’2  (“l  “in’ *3 

m  1 

-  (5-^r-lk-lW4,>2x7(“i-V],4 
C1  r 


(7.2.8) 


Subscript  1  indicates  contribution  from  \VY  •  Vri 

A. 

Subscript  2  indicates  contribution  from  C*Vy 

3C 


Subscript  3  indicates  contribution  from  to 


3p 


Subscript  4  indicates  contribution  from  |k  •  (•^xVoo) 

to  and  to.  are  the  vertical  motion  at  level  p  and  p. 
ml  ^m  ^1 

respectively . 

The  non-geostrophic  terms  will  be  approximated  by  the 
divergence  computed  from  the  geostrophic  part  of  the  equations 
for  the  preceeding  time  step. 

We  can  now  express  the  forecasting  equations  in  the 
following  formal  way  (compare  (2.2,  2.13  and  2.14). 

V2  {4+r  (^  +c,  )  }  =  F 

3t  rm  12  2Y1  mG 


+  F  +F  +F  +F 
ml  m2  m3  m4 


(7.2.9) 
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,2  8*1  a  f2  !!l  +  a  f2  !?2 

v  at  if  at  2f  at 


=  F 


1G 


+  F11  +  F12  +  F13  +  F14 

y2  9^2  .  .2  .  .  -2 

v  at  b2f  at  bif  at 


(7.2.10) 


F  +F  +F  +F  +F. 
2G  21  22  23  24 


(7.2.11) 


FmG-*  f1g  and  F2q  are  the  geostrophic  and  non-adiabatic 
terms.  They  correspond  to  the  right  hand  part  of  the  equa¬ 
tions  2.2,  2.13  and  2.14. 

The  non-geostrophic  terms  are  the  same  as  in  the  equa¬ 
tions  (7.2.6),  (7.2.7)  and  (7.2.8). 
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8.  NUMERICAL  SOLUTION  OF  THE  3-PARAMETER  MODEL 


8 . 1  GENERAL  ASPECTS 

The  equation  of  the  model  will  be  applied  on  a  polar- 

stereographic  projection.  The  polar-stereographic  plane  is 

assumed  to  cut  the  sphere  at  a  latitude  d>  .  In  the  numerical 

o 

computations  cf>  is  put  equal  to  60 °N,  however,  other  values 
of  can  easily  be  chosen.  The  grid-distance  d  can  also  be 
chosen  arbitrarily.  For  further  details  see  the  program 
description  in  Appendix  B. 

The  computational  area  can  have  any  form,  from  an  irregu¬ 
lar  octagon  to  a  square.  The  only  geometrical  condition  is 
that  the  inner  angles  of  the  area  must  be  90°  or  135°. 

The  coordinate  axis  of  the  grid  is  positively  oriented 
(see  Figure  3) .  The  computational  area  is  specified  by  the 

coordinates  of  the  corner  points  (x,/y  •••  x  /y  )  and  by  the 

-LX  o  o 

coordinates  of  the  north  pole  (x  ) . 


12 

li 


10 


9 


8 


7 


6 


xD/yn  ® 


p  jp 


5 


1 


2 


4 


3 


d 


d  1  2  3  4  5  6  7  8  9  10  11  12  13  14 


Figure  3.  Example  computational  area. 
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If  the  area  is  reduced  to  a  rectangle  x  =  x  ,  y  = 

1  2  1 


X3  =  X4'  y3  =  y4 '  etc* 


8.2  FINITE-DIFFERENCES 

In  order  to  transform  the  differential  equations  to 
finite-difference  equations  we  will  introduce  the  following 
finite-difference  notation  and  operators  (a  and  8  are 
arbitrary  quantities) : 


x  ^  iAs 


y  -  jAs 


t  -  t  At 


f(x,y,t)  fX.  . 

1  r  J 


(8.2.1a) 


J  (a  ,6)  -  — —  j)  (a  ,  8) 
4d 


(8.2.1b) 


3  a  A  Ta 
3t  eAt 


(8.2.1c) 


(8.2.2b) 


j,  o 

x=  0  yields  a2  -a  and  e=k 
x  =  h  yields  -  a°  and  e=l 
T  >  1  yields  at+1  -  at_1  and  e=2 


(8.2.2,g) 


34 


l+sincf) 


Introducing  the  map-scale  factor  m  =  — — —  n ^  and  inserting 


in  2 

the  quantity  y  =  (-5-)  we  get: 


V2(AiP1) 

fz 

-  al— 

(Aip2)  +  a2- 

fz 

y 

(i*2>  =  F1G, 

(8.2.3) 

V2  (A^2) 

-  b2  f2 

2  y 

(  At^2  )  +  bx- 

f2 

y 

(i+l>  -  F2G, 

(8.2.4) 

v2(AV 

-  %(Ail j  ) 
y  rM 

=  F 

mG. 

(8.2.5) 

q  is  an  empirical  constant  to  adjust  for  the  very  long  waves, 
-12  -2 

q  =  0.75*10  m  . 

Here  we  have  (for  simplicity  we  from  now  on  put  J|  =  J)  : 

F1G  =  V  +  V' 


2G 


F2 '  +  F2 " ' 


F  ' 


-  eAt*-(a3ws+a1H)  , 


T7I  II 
F1 


-  eAt^-  [J1+J2+f 2  (a2J4~a1J3)  ]  , 


F  1 

*2 


V 


-  eAti-  (b^+bjH)  , 

-  eAtj  [J5+J6+f2(b1J3-b2J4) ] , 


mG 

We  further  have : 


32  nm  2?lf 
^6  =  nm+2^2 ' 


eAt4  [J7+Jg+J9  4-jjCgWs  ]  . 


3-7  =  C-T)  —  C  „  £  +  C  „  , 

7  3  m  2  1  4  2  7  ' 

P8  =  -c2nm+c5?l- 
3 g  =  c4rim+CgC2+2c7f  , 
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J1 = 

J(V‘?i)# 

J2  = 

J(^-l;32)  , 

J3  = 

' 

il 

J(  W  , 

II 

LD 

l"0 

J(^m;C2)  ' 

J6  = 

J(^2;p6)  ' 

II 

r- 

l"0 

J(Ve7>  , 

00 

II 

JM;V ' 

C-I 

il 

J(^2;39)  , 

J10 

=  J(^m-2^1;ps)  , 

0)  = 
s 

4  J10_Cf  M"2^* 

Cf  is  an  empirical  constant  and  a  function  of  the  exchange 
coefficient  of  momentum  in  the  boundary  layer. 

8.3  COMPUTATION  OF  SENSIBLE  HEAT 

Equation  (4.1.4)  reads  in  finite-difference  form: 

hsens  =  °-5-io~2  h  <o.n^-((*xv 2+<v°)2>  +1)(Vh2 

over  ocean,  if  (T^-h^^-^)  >0  /  (8.3.1) 

H  =  0  over  land  and  if  (T  — h „ ip n )  <0  • 

DhjJNlb  S  /L 

Mo  =  [^(i+l)-2^(i+l)]-I^(i-l)-2^(i-l)] 

4y^o  =  j  +  D  “2^x(  j+D  ]~  j-D -2^1(j-D  ] 
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8.4  COMPUTATION  OF  LATENT  HEAT  (See  Gambo  1963) 

The  latent  heat  is  introduced  in  the  model  by  the 
expression : 


HLAT  -  -hL“*F* 

if 

“l  «-5l 

hlat  =  0 

if 

"l  >_51 

(8.4.1) 

hlat  -  0 

if 

eAt  tolerance 

0)’ 


-(61+62)  if  w1  <-(61+62) 


0). 


(8.4.2) 


U)x  = 


61+  62 


if  -  ( 6^+62)><:a)^-6 ^ 


6^  and  6^  are  here  two  tolerances  with  the  same  dimension 
as  r  is  the  precipitation,  and  <$2  is  introduced  for 

operational  purposes. 


E(T) 

£  L | 

F*  =  F* (p , T)  =  P 

c 

Lp 

pT2+|-^-E(T) 

P 


(8.4.3) 


E  (T) 


eL  f  1 
R 

o 


e  =  0.622 
L  =  2.5-106 


C  =  1004 
P 

T  =  273 
o 

E  =  0.611 
o 


P  -  P 

mean 


P  +P 
m  o 

2 


2f 

T  =  p  =  h6^i 

Rln(p^) 

m 
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8.4.1  Computation  of  V  and  D 


We  will  compute  the  humidity  in  the  layer 
Po  =  100  cb  Pt  =  30  cb  '  defining  as 
100 

V  =  7  O'  /  NV(P)  E(P)  dP 
30 

and  inserting  the  wind  profile  yields 

XV  =  emVm  +  elVVl  +  e2\V2  (8.4.4) 

where 

em~70’  a-^QQ+2  5E(70)  a^Q+20E(50)  a^Q+10E(30)  a.^  ]  , 

el=75'tl5E(100)b100+25E(70)b70  +  20E(50)b50+10E(30)b30]  ' 

e2=  W[15E(100)c100+25E(70)c70+20E(50)c50  +  10E(30)c30-1  * 

(8.4.5) 

The  constants  a^,  bp  and  c^,  where  p  assumes  the  values 
100,  70,  50,  30  are  computed  for  the  three  following 
alternatives : 


p  -$p^<p 


a  =  1 

P 

b  =  - 

P 


c  =0 
P 


2 (P-Pm) 


m 


P  ~P 
^o  rm 


II 

ppp<pm 


aP  =  1 
bp  =  ° 


c  = 


2  <  v-p’ 

Pm  Pi 


III 

0^P<P 


1 

p  “  p- 


a  =  & 


b  =  0 
P 


=  2p 


E ( 10 0 ) =2 . 5415 ,  E(70) =0 .9683,  E ( 50 ) =0 . 3527 ,  E(30)=0.0613 
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Correspondingly  we  get 


D  =  e  D  + 
m  m 


elDl 


+ 


S2D2 


(el+emC2)Dl  +  (e2-emcl)D; 


emc0w 
m  8  s 


8.4.2  Computation  of  Precipitation 

Equation  (4.2.6)  reads  in  finite-difference  form  (the  term 
VxVw  disregarded  and  w  replaced  with  q) : 

T+l  T-l  A  x 

q  =  q  +  eAt  H  , 


Ht  = 


^J('^T;qt)  ~  qT(DT 


+dw  T)  + 
s 


dif  f  ^ 


i2  T 

f  q 


(8.4.6) 


The  precipitation  is  computed  in  the  following  way  (precipita 
tion  is  indicated  by  r)  : 


qSAT  q(Ti|))  ~  q(h3^1+h4^2) 

If  (q  1-0.8  qSAT)  >0  then  qT+1  =  0.8  q 


T  4-  ] 

r  =  Ap(q  -0.8qSAT)  . 


T+l 


If  (q  0.8  qSAT)  <0  then  r  =  0  and  if 


(q  +1  0.2qSAT)  <0  then  q^  = 

r  is  accumulated  for  every  timestep  and 
prescribed  times.  See  Appendix  B. 


0.2 


qSAT  ‘ 


(8.4.7) 


printed  out  at  certain 
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8.5  NUMERICAL  SOLUTION  OF  THE  NON-GEOSTROPHIC  TERMS 


The  finite  difference  analogues  for  (7.2.9),  (7.2.10)  and 

(7.2.11)  read : 


V2  (AipJ  -  —  ( Ax/;  )  =  F  ^+—  Z  F  . 
rM  y  yM  mG  y  ^  mi 

2  2  4 

T2  (A^)  -  Bli  (A^)  +  a2i  (A#2)  =  F1G+MtjiFi. 


(8.5.1) 


2  2  4 

t  (A*2)  -  b2£  (A*2)  +  bl  Wy  =  F2G+£M  iIiP2 . 


Fml  2  ^c3  ^m+f  ^  c2^1+c4?2+c7f  ^ 


v_ 


a. 


(8.5.2) 


(8.5.3) 


+  ▼x1*V[-c2(Cm+f)+c5S1] 

V - ' 


a 


8 


+  Vx2-t[c4(Cni+f)+c6?2+2c7f]} 


(8.5.4a) 


Fm2  =  tDl(kl?l+k2?2+k3V  +  D2(k4?l+k5^2+k6V 


+  “s(k7?l+k8?2+k9?m)] 


(8.5.4b) 


Fm3  [Dl(k10?:l+kll?2+k12i;m)  +  °2  (k13?l+k14^2+k15?m) 


+  ws (k16Fl+k17?2+k18?m) ] 


(8.5.4c) 


Fm4  2^1  [k19Dl+k20D2+k21Ws] 


Y- 
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+  v^2*V'-k22Dl+k23D2+k24U3s-1 


v_ 


Y 


8 


+  ?i|>  •  V  [kocD, +k0,D0+k__a)  1 } 
rm  251  262  27  s 

V _ / 


Fn  =  - 

F12  -  IDl(k28!:r?m>+D2k29?l+"Sk30?l1 


F13  ^1 <k31Dl+k32D2+k33ws)  ' 


F14  “  !{V'<'l-T[k31Dl+k32D2+k33“sI} 


Y 


10 


21 


"  l{VXm-^2+Vx2'V(Sm+f+2C2)} 


F22  {Dlk34C2+D2  [k35?2  +  Wsk30C2* 


F23  ^2  (k36Dl+k37D2+k38a)s)  ^ 


?24  2^V^2,'i?*-k36Dl+k37D2+k38a)s-' 

v _ / 


Y 


11 


The  computation  of  the  non-geostrophic  forcing 

4  4  4 

2  F  .  *  E  F  .  and  ^AE  j  p  . 

*  i-1  ml  "  i-1  11  V  i=1  2i 


are  performed  by  a  separate  program.  See  Appendix  I 


(8.5.4d) 

(8.5.5a) 

(8.5.5b) 

(8.5.5c) 

(  8 . 5 . 5d) 

(8.5.6a) 

(8.5.6b) 

(8.5.6c) 

(8.5. 6d) 

function 
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8.5.1  Numerical  Coefficients  for  the  Non-geos trophic  Terms 
Expressions  for  constants  of  the  non— geostrophic  terms 

read : 


00  ^  = 
3 

t6Ws+t7Dl+t8D2 

00  = 
m 

t9Ws+t10Dl+tllD2 

“i = 

t12C0s+t13Dl+t14D2 

fc6  = 

k7  = 

z  6 

fc8  - 

Pi 

=  (c  -2)  — 

6 

k9 

:  e8(po',,.n) 

t10 

=  (c2-l) <p0-pm) 

kll 

=  -  c,  (p  -p  ) 

1  ^o  ^m 

t12 

=  c  !i 

8  2 

^13 

Pi 

=  “C2  ~ 

2 

k14 

Pi 

=  (c  -2)  — i- 
2 

kl  = 

C2  ^C2~3^ 

k2  - 

"C2C4 

k3  = 

C2C7 
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k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 


4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


~C1C2 


Cn  C  .  ” c ^ 

14  6 


c?  ( 2— c-j^) 


“C2C8 


C4C8 


C8(1~C7) 


2t2C8 

-2 c 

_t7C8 

“2t3C8 

2(t5"t8)c8 


t8C8 


2tlC8 


2(t4  tg)c8 


2t2C8 

_2t3C8 


8 
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k21  “  2tlC8 


k22  ~  “  2  (t3_ht 7)  c8 


k23  "  2(t5  t8)°8 


k24  "  2 (t4  t6)  C8 


k25  t7C8 


k26  -t8°8 


k27  t6C8 


k  =  2-c„ 
K28  2 


V  ss  4.  /^ 

K29  "1 


k30  =  +C 


10 


"31 


‘32 


‘33 


(PcTlV 


'11 


(po 

(Po'pm) 


k34  =  ~c; 


‘35 


c^-2 


"36 


k10  k13 

(Pm"pl) 
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1r 

37 

kll  fc14 

"  'fm-Pl1 

k38 

t9 t12 

For  po=100,  Pm=50  and  P1=30  we  have  the  following 
numerical  values  for  the  constants : 


fc6 

0.0588235 

tl 

-2.941175 

-7.058825 

fc9  = 

0.411765 

t10 

=  -20.588250 

i — 1 
i — 1 

•P 

=  -29.411750 

k12 

=  0.176471 

k13 

=  -8.823525 

fc14 

=  -21.176475 

II 

i — 1 

-0.438291 

k2  = 

-0.276816 

k3  = 

0.034602 

k4  = 

-0.346020 

k5  = 

-0.507498 

k6  = 

0.083045 

11 

-0.006920 

k8  = 

0.005536 

5* 

11 

0.011073 

kio 

=  -0.438294 

kll 

=  -0.276817 

45 


k12 

0.034602 

k13 

- 

-0.346021 

k14 

- 

-0.507498 

k15 

= 

0.083045 

k16 

= 

0.016609 

k17 

= 

0.005536 

k18 

= 

-0.000692 

k19 

= 

-0.438293 

k20 

= 

-0.346021 

k21 

= 

0.016609 

k22 

= 

-0.276817 

k23 

= 

-0.507497 

k24 

= 

0.005536 

k25 

= 

0.034602 

to 

CTs 

= 

0.083045 

k27 

- 

-0.000692 

k28 

= 

] .411765 

k29 

= 

+0.588235 

k30 

SB 

+0.0117647 

k31 

= 

0.411765 

k32 

= 

0.588235 

k33 

= 

0.0117647 

k34 

= 

-0.588235 

k35 

= 

-1.411765 
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-0.588235 


k3?  =  -0.411765 
k3g  =  0.0117647 


47 


9.  INITIALIZATION 


The  stream  functions  \p  are  computed  from  the  geopotential 
Z  by  the  relation 

gV2Z  =  V-  (fVi|>)  or 
V2\p  =  |V2Z  -  i-Vf*Vi|); 
inserting 

Vip  =  ^VZ  yields 

V2\p  =  |v2Z  -  2_Vf-VZ. 

The  boundary  values  of  ^  are  computed  by  the  relation  (y  is  a 
constant) 


(9.1) 


(9.2) 


3^ 

3s 


g  3_Z 
f  3s 


+  T 


(9.3) 


<  i  . ,  ,  1  g  3Z, 

easily  seen  that  y  =  -—  gs  — ds 


If  we  assume  that  there  is  no  net  flow  out  of  the  area  it  is 

The  integration  is  performed 
along  the  boundary  of  the  area  in  positive  order.  The  integra¬ 
tion  starts  in  point  x^/y^  where  we  put 


^(x1;y1)  = 


fuTTTT)  z(xi;yi) 


(9.4) 


Z  is  computed  from  \p  in  a  similar  way  using  the  following 
equation: 

V2z  =  -v2^  +  -Vf«Vi|>.  (9.5) 

9  9 


Boundary  values  for  Z  are  computed  initially  and  stored  in  a 
special  string.  See  Appendix  B. 
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9.1  NUMERICAL  SOLUTION 


Equation  (9.2)  r^ads  in  finite  difference  form: 

V2\p  =  |-v2z - Vf •  VZ 

1  2  f  2 


(9.6) 


where  the  standard  five  point  formulas  are  used  to  compute 
V2Z  and  VfVZ  is  defined  as 

Vffz  +  (fi-Vi’j'VVi’j 


Equation  (9.3)  reads  in  finite  difference  form: 

\+l 


=  ^  +  2gfert  +  Y(is>k 


where 


i  ,a  zk+i~zk 

T  =  E  9  k-i  Wfk 


(9.7) 


(9.8) 


Figure  4  defines  k  and  the  order  of  integration. 


Z  to  ^  and  to 
Z  is  computed  by 
the  program  STREAMF 
(see  Appendix  B) . 


Figure  4.  Definition  of  k  and  order  of  integration 
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9.2  INITIALIZATION  OF  THE  SPECIFIC  HUMIDITY 

Since  we  are  using  the  ^-functions  as  the  time-dependent 
^^^i^hles f  we  have  to  modify  the  initial  humidities  in  order 
to  make  them  consistent  with  the  time-integration. 

Therefore,  we  put 


(9.2.1) 


n>  _  2g  ,Z500  Z1000,_g 
Z  Rln2  v  2  >~£ 


^3Z1+^4Z2^  '  an<^ 


o 


n>  =  hsW2* 
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10.  TIME-INTEGRATION  AND  TREATMENT  OF  LATERAL  BOUNDARY  VALUES 


The  forecast  is  basically  computed  in  6-hour  intervals 
but  this  interval  can  easily  be  changed.  The  first  time  step 
in  each  interval  is  non-centered .  Smoothing,  elliptization 
and  printing  of  results  (if  so  desired)  are  performed  at  the 
end  of  every  interval. 


N 


a^ 

at 

ND 

kT 

N 


Ai p_ 

eAt 


At  =  1  hour  in  this  example 


number  of  At  for  the  interval 

time  step  index;  k  =  1,2,3...ND+1 

index  for  every  interval  integration  (e.g. 
interval) 


(10.1) 


6 -hour 
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Initial  height  fields  Z°  are  stored  on  secondary  storage 
during  the  whole  computation.  In  the  case  of  variable 
boundary  conditions,  Z°  is  followed  by  ZN  (N=l,2,3...) 
(forecasts  for  each  interval) . 

Assuming 


lt=Tlf  or  (ZN-Z°) 


(10.2) 


the  stream  function  at  time  step  k  can  be  interpolated  from: 


1 Pk  =  (1>°  ~  I  Z°)  +  I  [  (1-a,  )  ZN  +  ol  ZN+1] 


(10.3) 


where 

“k 


kT-1 

ND 


At  each  time  step  this  interpolated  stream  function  is  mixed 

k 

with  the  forecasted  stream  function  ip  in  the  following 

prog  3 


way : 


Boundary  values 

1  grid  point  inside 
the  boundary 

2  grid  points  inside 
the  boundary 

3  grid  points  or  more 


ij>  .  =  ip^  +  wl(fk-i|)k  ) 
mod  rprog  r  rprog 


\p  .  =  \p^  +  w2  (ip^-rp^  ) 

ymod  rprog  rprog 


ip  ,  =  ip^  +  w3(^^-iJjK  ) 
mod  prog  prog 


inside  the  boundary  ip  .  =  ip  , .  „  . . 

J  rmod  rprog  (10.4) 

wl,  w2  and  w3  are  3  predetermined  constants  with  o=Swl, 
w2 ,  w3^1 

The  following  flow  diagram  illustrates  the  treatment  of 
the  boundary  values: 
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Flow  diagram:  Treatment  of  boundary  values 


/" 
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II.  THE  CRITERION  OF  ELLIPTICITV 


An  iterative  procedure  is  used  to  modify  a  stream- 
function  field  so  that  the  ellipticity  criterion 

V2ip  +  |  >  0 

is  valid  in  all  points  of  the  field. 

Each  point  is  tested  with  the  following  formula: 

+  |  -  e  =  6 

where 


2 

▼  '  ip 


=  v 


i+l,  j 


*i- 


1/ j 


+ 


*i,j+i  + 


and 


e  =  0.001-f. 

If  6  <  0,  ip.  .is  modified  by 

^i ,  j  “  ^i ,  j  -  2^jT  where  k  is  a  convergence  parameter. 

This  means  that  the  vorticity  increases  by  2(e-S)*k  in 
the  point  (i,j)  and  decreases  by  0.5  (e-6) *k  in  the  four 
surrounding  points  (i+l,j),  (i-l,j),  (i,j+l)  and  (i,j-l). 

This  procedure  guarantees  that  6  becomes  positive  in  the  point 
(if j)  r  but  not  necessarily  in  the  surrounding  points.  The 
test  must  therefore  be  repeated  for  all  points  until  the 
criterion  is  valid  in  the  whole  field.  With  a  suitable  choice 
of  k  the  method  is  convergent. 

k  can  be  found  empirically  and  is  estimated  to  be  of  the 
order  k  -  0.85. 
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APPENDIXES 
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APPENDIX  A 


GENERAL  PRINCIPLES 
IN  THE 

PROGRAMMING  OF  THE  THREE  PARAMETER  MODEL 

The  programming  has  been  performed  in  a  very  general 
way.  It  is  therefore  possible: 

a.  to  generate  initial  fields  and  run  the  model  in  a 
channel  with  variable  cyclic  boundary  conditions  (this  is  the 
case  for  the  actual  program) . 

b.  to  run  the  model  with  constant  boundary  values  on 
a  given  area  (see  above) . 

c.  to  run  the  model  with  variable  boundary  values,  the 
boundary  values  taken  from  another  model. 

The  model  needs  13  fields  in  the  core ; 


1  field  for  the  coriolis  parameter  F, 

2 

1  field  for  y  =  MY, 

d 

1  special  field  indicator  MARK, 

10  fields  for  the  model  computations  F1-F10. 


In  addition  to  this,  40  fields  (28  for  the  quasi-geostrophic 
part  of  the  model)  are  specified  on  secondary  memories  (disks, 
drums  or  large  extended  core) . 

"Household"  programs : 

PUT1:  This  program  computes  FPAR  (see  common, 

Subroutine  JMIMA) 
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MARKF : 

This  program  computes  a  special  integer  field 

MARK.  The  program  MARK  specifies  status  of 

the  field.  Points  outside  the  area  =  0,  points 

inside  the  are  <0  and  points  on  the  boundary 

>0.  The  corner  points  etc.  are  specified 

according  to  given  examples. 

MYFF  : 

This  program  computes  f  and  y  and  puts  the 

result  in  the  fields  F  and  MY. 

RANWT : 

Writes  fields  on  secondary  storage  e.g.,  disk 

drums,  or  extended  core  storage. 

RANRD: 

Reads  fields  from  secondary  storage  e.g., 

disk,  drums,  or  extended  core  storage. 

MAP  : 

Prints  pattern  on  line  printer; 0  (zero)  points, 

resolution  and  map  scale  are  specified. 

GENCH : 

Generates  initial  state  for  a  channel  flow. 

BMOVE : 

Administrates  computation  of  cyclic  boundary 

conditions  for  a  channel  flow- 

Since  the  axis  of  the  grid  is  defined  differently  from 
what  is  used  in  Fortran,  the  programming  of  finite  difference 
operators  should  be  performed  in  a  special  way.  This  is  not 
necessary,  but  speeds  up  the  computation  considerably.  As 
an  example,  subroutine  JACOB  reads: 
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MI  =  M-l 
D<j>  10  I  =  2, MI 
J1  =  JMIN (I) +1 
J2  =  JMAX (I) -1 
K  =  (J1-1)*M+I 
D<f>  10  J  =  J1,J2 

C  (K)  =  (A (K+l)  -  A(K-l)  *  (B  (K+M)  -  B  (K-M)  ) 
10  K  =  K  +  M 
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APPENDIX  B 


PROGRAM  SPECIFICATIONS 

The  following  programs  are  written  for  the  3-parameter 
model. 


Level 

1: 

Main  program 

Level 

2  : 

Subroutines 

called  by  level 

Level 

3: 

Subroutines 

called  by  level 

Level 

4  : 

Subroutines 

called  by  level 

Level 

Program 

1 

PR0G3P 

2 

PUT1 

3 

JMIMA 

2 

C0EFF3P 

2 

MARKF 

2 

MYFF 

2 

STREAMF 

3 

BDRGDR 

3 

BDRVAL 

2 

SATUR 

2 

RANWT  ( RANRD) 

2 

ELLIPT 

2 

GENCH 

2 

ASMUT 

2 

MAP3P 

3 

MAP 

2 

STEP3P 

61 


Level 

3 

3 

3 

3 

3 

3 

3 

4 
4 


Program 

JACOB 

ABSVOR  (RELVOR) 
MIXF 

HELM  (POIS) 

HELMSYS 

BMOVE 

STEPEXT 

GRADPR 

VELPOT 


62 


Program  PR0G3  (Main  program  for  3-parameter  model,  see  flow 
.  ■  ■—  diagram  B-l  and  description.) 

Arrays 

Name  Dimensions  Description 

F1-F10  2  Working  fields  in  core.  See  flow 

diagram  B-l. 

F ,MY  (real)  2  Fields  for  Coriolis  parameter,  f, 

and  y  =  (2.)  z .  m  is  the  mapscale 
factor  for  a  polar  stereographic 
projection  and  d  is  the  grid 
length.  F  and  MY  are  generated 
by  the  subroutine  MYFF. 

MARK  2  Marking  of  each  grid  point  by  a 

special  integer  (index) .  See 
description.  MARK  is  generated 
by  the  subroutine  MARKF. 


UPS,  UPM,  UP1  1  Zonal  wind  profiles  at  the  levels 

p  ,  p  and  p.  for  the  initial 
fields  in  the  channel  case.  The 
values  are  introduced  by  DATA 
statements . 

WX  1  Working  array  for  generation  of 

initial  fields  in  the  channel 
case.  Dimension  shall  be  at 
least  equal  to  the  number  of  grid- 
points  across  the  channel. 
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Name 


Dimensions 


Description 


NX,  NY 


1 


PSIC,  PSIS,  1 

LAMC,  LAMS 


ZB , PSIB ,FB , SB  1 


IT 


2 


Wave  numbers  for  perturbations 
of  the  initial  fields  in  the 
channel  case  in  x  and  y  direc¬ 
tions.  The  values  are  introduced 
by  DATA  statements. 

Amplitudes  (PSI)  and  phase  lags 
(LAM)  perturbations  of  the  initial 
field  with  cosine  and  sine  profiles 
in  the  y  direction  respectively. 

For  both  profiles  a  sine  wave  is 
used  in  the  x-direction.  The 
values  are  introduced  by  DATA 
statements.  See  further  descrip¬ 
tion  of  GENCH. 

Boundary  strings  for  storage  of 
boundary  values  in  counterclock¬ 
wise  order,  starting  with  corner 
point  1  (see  FPAR) .  Used  by  the 
subroutine  STREAMF . 

Storage  of  the  number  of  iterations 
for  the  solution  of  a) :  the  system 
of  two  Helmholz  equations  and  b) : 
the  Helmholz  equation  for  the 
mean  field.  Values  for  an  inte¬ 
gration  interval  are  stored  in  the 
array  by  the  subroutine  STEP3P . 
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Name 


Dimensions 


MAPHOUR  (1)  1 

NSMUTT  (1)  1 

MELLIPT  (1)  1 

LETACC  (1)  1 


Description 

Integers  giving  the  hours  for 
mapping  (2) . 

Integers  giving  the  hours  for 
smoothing  (2) . 

Integers  giving  the  hours  for 
ellipticity  (2) . 

Integers  giving  the  hours  at 
which  the  accumulation  of 
precipitation  shall  be 
interrupted  (2) . 


^  The  values  shall  be  given  in  DATA  statements  and  must  be 
multiples  of  the  integration  interval  (e.g.,  6  hours).  Unused 
positions  are  indicated  by  -1. 

( 2) 

This  value  can  be  increased  if  this  is  necessary. 
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Flow  Diagram  B-l:  PR0G3P 
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Common  areas 


/FPAR/ 

Name 

TYPe 

Description 

IC  (  8)  ,  JC  ( 8) 

integer 

i-  and 

j -coordinates 

for  the 

corner 

points  of  the 

area. 

XPOL , YPOL 

real 

i-  and 

j -coordinates 

for  the 

north  pole.  Specified  relative 
to  the  origin  of  the  grid. 

R  real  radius  of  the  earth. 

RE  real  Distance  from  the  pole  to  the 

equator  on  the  map  in  grid-length 
units . 


DS 

real 

grid- length  (in  meters) . 

JMIN(IOO) , 
JMAX(IOO) 

integer 

minimum  and  maximum  j -coordinate 

for  each  i-column  in  the  field. 

M,N 

integer 

Dimension  of  the  field  arrays. 

KIND 

integer 

Channel  indicator.  Channel  =  1, 

No  channel  =  0.  Value  shall  be 

given  by  DATA  statement. 

/F0RM1/ 

IC1  (8)  , JC1(8) 

integer 

Same  as  IC(8),JC(8)  for  actual 

field.  Values  introduced  by  DATA 
statement . 
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Name 


JMINl(lOO)  , 
JMAXl(lOO) 

XP,YP,D1 


/KANAL/ 

FIM 


/DRM/* 

MN 

NDIM 


NFLD 

SL 


Type  Description 

integer  Same  as  JMIN(IOO),  JMAX(IOO) 

real  Same  as  XPOL,  YPOL,  DS .  Values  shall 

be  introduced  by  DATA  statements. 

Dl  is  given  in  km. 


real  Used  for  channel  computations . 

Latitude  for  the  middle  of  the  channel. 
Used  for  computation  of  B-plane.  To 
be  given  by  DATA  statement. 


integer  =M*N.  To  be  computed  in  main  program 
integer  Core  memory  space  available  for  field 
arrays.  The  arrays  F1-F10 ,  F, MY, MARK 
must  be  included  in  this  area. 
Additional  space  is  used  for  storage 
of  fields  in  COMMON  area//, 
integer  Number  of  field  arrays  in  core 
real  Working  area  for  core  memory  fields . 

The  dimension  must  be  at  least  SL(NDIM) 


Arrays  stored  at  extended  core,  disks  or  drums.  1  in  the  end 
of  the  name  indicate  timestep  x,  2  in  the  end  of  the  name 
timestep  x-1. 


* 

Not  necessary  to  specify  if  only  extended  core  storage  is  used. 
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PSIM1 ,  PSI11 , PS I 21 
PSIM2 , PS I 12 ,PSI22 


HUM1 ,HUM2 , DIV1 
DIV2 ,WS ,HEAT 

J789  ,J12,J56,J3 
PS ,TS ,PREC 


STRM , STR1 , STR2 , ZM1 
Zll , Z21 , ZM2 ,  Z12 ,  Z22 

H11,H21,HM1 


H12 ,H22 ,HM2 


H13 ,H2  3 ,HM3 


VI ,V2 ,VM 


ID (200 ) * 


Not  necessary 
is  used. 


real 


real 


Storage  arrays  for  stream  func¬ 
tions.  See  flow  diagram  B-l. 
Storage  arrays  for  humidity, 
divergence,  vertical  velocity 


at  the  lower  boundary,  and  heat* 
Storage  arrays  for  Jacobians. 
Storage  arrays  for  surface 
standard  pressure,  sea  surface 
temperature,  and  accumulated 
precipitation. 

Storage  arrays.  See  Chapter  10 
(lateral  boundary  mixing). 

Fields  for  advection  of  vorticity 

by  the  divergent  wind,  thickness 

field  (1,2)  and  mean  field  (M) 

3 

W¥t+  ?V*\V  for  thickness  fields 
(1,2)  and  mean  field  (M) . 

Twisting  term  for  thickness  fields 
(1,2)  and  for  mean  field,  M. 
Velocity  potential  for  thickness 
fields  (1,2)  and  for  mean  field 


(M)  . 


integer  Catalog  array  for  direct  memory 
access . 


to  specify  if  only  extended  core  storage 
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/COEFF/ 

A1,A2 . ,T5 


/COEFF2/ 

T6  ,  TF . ,T38 


/RUNPAR/ 

DELT 


NTSTEP 


ALFASYS,  ALFAM 
ALFAZ , ALFAPSI , 
RESSYS , RESM, 
RESZ , RESPSI 


real  Constants  used  in  the  model. 

Computed  by  the  subroutine 
C0EFF3P  from  given  values  on 
stability  and  pressure  levels. 

real  Constants  used  for  the  computation 

of  the  non-geostrophic  terms. 
Computed  by  C0EFF3P  from  given 
values  on  stability  and  pressure 
levels . 


real  Timestep  in  sec.  computed  in  the 

main  program. 

integer  Number  of  timesteps  for  an  inte¬ 
gration  interval  (6  hours) . 

To  be  defined  in  a  DATA  statement, 
real  Overrelaxation  coefficients  (ALFA) 

and  maximum  residual  in  the  solu¬ 
tion  for  the  system  of  the  two 
Helmholz  equations  (SYS) ,  Helmholz 
equation  for  the  mean  field  equa¬ 
tion  (M)  ,  solution  of  Z  from  \p 
(Z)  and  solution  of  ^  from  Z 
(PS I) .  To  be  defined  by  a  DATA 
statement. 


70 


Q , FOCEAN , FCONT 

real 

The  Helmholz  term  for  the  mean 

field  equation,  friction  coeffi¬ 
cients  over  ocean  and  land.  To 

be  defined  by  a  DATA  statement. 

WGT1 ,WGT2 ,WGT3 

real 

Weights  for  mixing  near  the 

boundary  of  boundary  fields  with 

forecasted  fields.  To  be  defined 

by  a  DATA  statement. 

ADIFF 

real 

Diffusion  coefficient  for  humidity; 

it  is  defined  by  a  DATA  statement 

Parameters  only 

specified  in 

Data  statements. 

STAB1 

real 

Static  stability  of  layer  1 

STAB2 

real 

Static  stability  of  layer  2 

PNIVS* 

real 

Pressure  level  PQ 

PNIVM 

real 

Pressure  level  Pm 

PNIV1 

real 

Pressure  level  Pi 

I  VAR 

integer 

Indicates  if  variable  lateral 

boundaries  are  used;  IVAR=0,  constant 

boundary  values;  IVAR=1,  variable 

boundary  values . 

KIND 

integer 

Indicates  if  channel  is  used 

(cyclic  boundary  conditions); 

KIND=1,  channel; 

KIND=0,  no  channel  (polar  stereo¬ 
graphic  projection)  . 
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Subroutine  STEP3P (F1,F2,F3,F4 ,F5 ,F6 ,F7 ,F8 ,F9 ,F10 , MY, F, MARK, 

IT , IVAR,M,N) 

The  subroutine  performs  a  computation  for  a  time  interval 
(6  hours)  by  the  3P-model  including  humidity  and  precipitation 
prediction.  (See  Flow  Diagram  B-2.) 


Subroutine  Parameters 

Description 

F1-F10 

Working  fields  . 

MY 

^-parameter  field  . 

F 

Coriolis  parameter  field  . 

MARK 

Field  indicator  field. 

IT 

Array  to  store  number  of 

iterations  for  every  timestep . 

IVAR 

Indicator  for  constant  or 

variable  boundary  conditions  ; 

* 

IVAR  =  0  constant  boundary  ; 

IVAR  =  1  variable  boundary. 

M,N 

Field  vector. 

Parameters  specified  only  in  DATA  statements 


RGAS 

R  =  287 

EE 

E  =  0.622 

HL 

L  =  2.5*106 

CP 

Cp  =  1004 

TO 

To  =  273 

EO 

Eo  =  0.611 

DELI 

61 

DEL2 

<S2 

TOL 

Tolerance 

For  explanation  see  chapter  8.4. 
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Subroutine  MARK 


The  subroutine  mixes  field  FA  (corresponding  to 
limited  area)  with  FB  (corresponding  to  a  field  taken 
from  a  large  area)  . 


large  comp, 
area 


FA 


limited 
comp,  area 


W1'W2'W3 

MN 


Weight  factors 
Field  vectors 
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MARK- fie Id 


The  MARK-field  generated  by  the  subroutine  MARKF 
can  be  described  by  the  following  examples. 


1 .  Channel  field. 

0  10  10  10  10  10  10  10  10  10  10  10  10  10  10  10  10  10  10  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -l  -l  -l  -l  -1  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -l  -1  -1  -l  -l  o 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -l  -l  -l  -1  -1  o 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  o 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -l  -1  -1  -1  -1  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -l  -l  -1  -1  -1  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  _l  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  0 

j  0-1-1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  0 

0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  _1  0 

it  0  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  0 

02222222222222222220 

- *»-  i 

2 .  Ordinary  octagonal  field. 
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3.  Fields  with  2  rectangular  points 


3 

1 1 


20 

10 

10 

10  10  10  10  10  10  10  10  10  10  o' 

14 

-9 

-8 

1 

00 

t 

CO 

1 

00 

I 

CO 

1 

00 

1 

00 

I 

CO 

1 

00 

1 

00 

1 

wj 

14 

-9 

-10- 

■lo-Io-io-io-io-io-io-io-io-io-id' 

14 

-9 

-10 

-1  -1  “1  -1  -1  -1  -1  -1  -1  -1  -t 

14 

-9 

-10 

-1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1 

14 

-9 

-10 

-1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1 

14 

-9 

-10 

-1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1 

14 

-9 

-10 

-1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1 

14 

-9 

-10 

-1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1 

14 

“9 

-10 

-1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1 

14 

-9 

-10 

-1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1 

15 

^10 

lei  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1 

1  -1  -1  -1  -1  -1  -1  -1  -1  -3 

i  -i  -i  -i  -i  -i  -i  -i  -i  -: 

-l  -l  -l  -l  -l  -i  -l  -l  -i 

0-10-10-10-10-10-10-10-10-10 

-6  -6  -6  -6  -6  -6  -6  -6  -6 


2  18 1 


3 

i  i 


0 

0 

0 

0 

<fl 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

19 

0 

0 

0 

<4 

-8 

-8 

-8 

-8 

-8 

-8 

-8 

-8 

-8 

-8 

-8 

-8 

-7 

6 

0 

0 

ao-io- 

-10-10- 

-10-10- 

-10- 

•10- 

-10- 

-10- 

-10- 

-10 

-10 

-7 

6 

0 

y* 

-i 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-10 

-7 

6 

<L 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-10 

-7 

6 

r!3 

\S 

jrio' 

l<L 

-1 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-10 

-7 

6 

14 

-9 

-10 

-1 

-1 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-10 

-7 

6 

14 

-9| 

-10 

-1 

-1 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

n 

-7 

6 

14 

-10 

-1 

-1 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-10 

-7 

6 

14 

:-9 1 

-10 

1-1 

-1 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-10 

-7 

6 

14 

-9 

-10 

-1 

-1 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-10 

-7 

6 

14 

-9 

-10 

-1 

-1 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1| 

h-10 

-7 

6 

14 

-9 

-10 

-1 

-1 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-lj 

-id 

-2J 

'  J 

14 

-9 

-10 

-1 

-1 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

14 

-9 

-10 

-1 

-1 

-l 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-3s 

'iC- 

0 

14 

-9 

r10- 

-10- 

-10- 

-10- 

-10- 

-10- 

-10- 

-10- 

-10- 

-10- 

•10- 

-10- 

-10,. 

0 

0 

14 

-9 

-6 

-6 

-6 

-6 

-6 

-6 

-6 

-6 

-6 

-6 

-6 

-6 

-2^ 

0 

0 

0 

17 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

3" 

0 

0 

0 

0 

i 


75 


Subroutine  RANWT  (A, ALFA) 


This  subroutine  writes  field  A  (in  core)  to  field  ALFA 
(in  secondary  storage). 

Subroutine  RANRD  (A, ALFA) 

This  subroutine  reads  field  ALFA  (in  secondary  storage) 
into  field  A  (in  core) . 

Subroutine  MIXF (FA, FB ,MARK ,W1 ,W2 ,W3 ,M,N) 

This  subroutine  mixes  field  FA  (corresponding  to  limited 
area)  with  FB  (corresponding  to  a  field  taken  from  a  large 
area) . 
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Subroutine  HELMSYS  ( Zl , Z2 ,F0RC1 ,F0RC2 ,F ,MY ,FMY , A1 ,A2 ,B1 , B2 , 
ALFA ,  RESIDUE  ,  IT  , ]VL,N) 

This  subroutine  solves  the  following  system  of  Helmholz 

equations  by  relaxation: 

2  2 

V2(Z1)  -  Al^  (Zl)  +  A2^ —  ( Z 2 )  =  F0RC1 
2  2 

V2(Z2)  -  B2^—  ( Z 2 )  +  Bl^ —  (Zl)  =  F0RC2 


Subroutine 

Parameter 

Zl ,  Z2 


F0RC1 ,F0RC2 
F 

MY 

FMY 

A1,A2,A3,A4 

ALPHA 

RESIDUE 

IT 

M,N 


Description 

Fields  to  be  solved  by  relaxation. 
First  guesses  in  Z1,Z2  before 
using  subroutine . 

Forcing  functions. 

Coriolis  parameter  field  . 
la-parameter  field. 

Working  field . 

Physical  parameters  (constants)  . 
a-overrelaxation  coefficient. 
Residual,  R,  in  the  solution 
of  the  system. 

Number  of  iterations . 

Field  vectors . 


77 


Subroutine  MAP3P (II , 12 , 13 , 14 , 15 ,16 , 17 , 18 ,19 , I1Q ,111 ,F1 , F2 , 
F3,F4 ,F5 ,F6 ,F7 , F8 ,F9 ,F10 ,F,MARK ,M,N,IDAY,ITIME,NTIME, PNIVS , 


PNIVM ,  PNIV1 ,  ZB1 ,  ZB2  ,  ZB3) 

Printing  on  line-printer  of  forecast  fields  in  zebra 
patterns.  Heights  are  computed  from  stream  functions. 


S  lib  routine  Parameters 

Description 

11-111 

Indicators.  I  =  0 :  no  printing 

1=1:  printing.  The  numbers 

refer  to  the  following  fields  : 

F1-F10 

XI:  Surface  pressure. 

12:  Height  for  level  p 

13:  Height  for  level  p™  . 

14:  Thickness  (p  -p  ).  "L 

15:  Lower  vertical  velocity  w..  . 

16 :  Precipi table  water  .  1 

17:  Accumulated  precipitation. 

18:  Relative  humidity. 

19:  Stream  function  for  level  p  . 

110 :  Stream  function  for  level  ps  . 

Ill:  Stream  function  for  level  p™  . 

Working  fields.  (See  flow  diagram  B-2) 

F 

Coriolis  parameter  field  . 

MARK 

Indicator  field  (see  MARKF)  . 

M,N 

Field  dimension  . 

I  DAY 

Year,  month  and  day  as  one  integer. 

ITIME ,NTIME 

Initial  and  forecast  time  . 

PNIVS , PNIVM, PNI VI 

Pressure  levels  in  the  model  . 

ZB1,Z.B2,£B3 

Working  arrays  for  boundary  strings  . 
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Flow  Diagram  B-2:  MAP3P 
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Subroutine  STREAMF  ( Z , PS I , R , F ,  MARK ,  ZB,PSIB,FB  , SB, GAMMA, IND, 


IT ,M,N) 


The  subroutine  computes 

the  linearized  balance  equation. 

Subroutine  Parameter 

Description 

z 

Z  field. 

PS  I 

PSI  field. 

R 

Forcing  function  for  the 

Poisson  equation  in  the 

2 

solution  of  V  \p  =  F(Z)  or 

V2Z  =  Z  (ip)  . 

F 

Coriolis  parameter. 

MARK 

MARK-field. 

ZB 

String  of  boundary  values  for  Z^ . 

PS  IB 

String  of  boundary  values  for  1^. 

FB 

String  of  boundary  values  for  f^_ 

SB 

String  of  boundary  values  for 

GAMMA 

<AS)k. 

Y- 

IND 

See  subroutine  comments. 

IT 

Number  of  iterations  for  solving 

the  Poisson  equation  by  relaxation 

M,N 

Field  vectors. 
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Subroutines  called  by  STREAMF 


BDRGRD  (SB) 

Computes  SB. 

BRDVAL (A,AB ,M,IND) 

Computes  boundary  values  from 

field  A  and  stores  the  boundary 

values  in  string  AB  or  reverse: 

IND  =0,  AB  =  A;  IND  =  1,  A  =  AB 

M 

Field  parameter 

IND 

See  subroutine  comments 

POIS  is  an  entry  point  in  HELM  for  the  solution  of  a  Helmholz 
equation  by  Liebmann  relaxation. 

HELM (Z ,FORC,Q, ALFA, RESIDUE, IT, M,N) 


Subroutine  Parameter 

Description 

Z 

Z-field. 

FORC 

Forcing  function. 

Q 

Helmholz  coefficient. 

ALFA 

Overrelaxation  coefficient. 

RESIDUE 

Residual . 

M,N 

Field  vectors . 

IT 

Counts  the  number  of  iterations 

required  to  obtain  a  solution. 
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Subroutine  SATUR(TM,QSAT) 


This  subroutine  computes  QSAT ,  integrated  mixing  ratio 

at  saturation  as  a  function  of  the  mean  temperature  between 

500  and  1000  mb.  TM  values  of  QSAT  are  given  for  each  whole 

degree  between  -50 °C  to  +20°C  in  the  subroutine.  Linear 

interpolation  between  whole  degrees  is  employed  and  the  result 

2 

QSAT  is  given  in  TON/m  /cb . 
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Subroutine  MAP (Q, M, N, PS ,QZ,QD, SCALE) 

Prints  a  field  in  "zebra  pattern"  on  line— printer . 
Subroutine  Parameters  Description 


Q 

M,N 

DS 

QZ 


QD 


Field  to  be  printed. 

Field  vectors. 

Grid  distance  (meters) 

Isoline  corresponding  to  000, the 
line  towards  999  on  printer 
(indicated  by  heavy  line  below) . 
Resolution  indicator  (see  below) . 


000000000  ) 
000000000  / 


Illllllll 

111111111 


QD 

QD 


SCALE 


222222222 

222222222 

Map  scale  factor, 
(unit  106  m) 
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Subroutine  JACOB (A, B ,C,M,N) 


This  subroutine  computes,  a  Jacobian  operator  (of  the  form) 


C  -  J(A,B)  (Ai+i  Ai-l*j(Bj+l  Bj-l)i 
(Aj+l~Aj-l}  i  CBi+l~Bi-l)  j 


Subroutine  parameters  Description 

A,B,C  Fields  according  to  formula. 

M,N  Field  vectors  . 
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Subroutine  ABSVOR(PSI ,VOR,F,MY ,MARK,M»N) 


This  subroutine  computes 
vorticity . 

a.  ri  =  yV2iJ>  +  f 

b.  C  =  yv2^ 

Subroutine  Parameters 
PS  I 
VOR 
F 

MY 

MARK 

M,N 


the  relative  or  the  absolute 

ABSVOR 

RELVOR  (via  entry  point) 

Description 

ip  field  . 
n  or  £  field  . 

Coriolis  parameter  field  . 
y  parameter  field  . 

Field  indi cator  array . 

Field  vectors  . 
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Subroutine  GENCHjPSI ,M,N,PSIO,U,NU,NWAVE,NX,NY,PSIC,LAMC, 
PS IS , LAMS ,WX) 


Subroutine  Parameters 

Description 

PSI 

Result  field. 

M,N 

Field  vectors. 

PSIO 

Constant  value  (^Q)  . 

U 

Zonal  wind  speed. 

NU 

Resolution  of  zonal  wind  speed. 

NWAVE 

Number  of  waves. 

NX 

Wave  numbers  as  a  function  of 

channel  length  in  x-direction  (nx)  . 

NY 

Wave  numbers  as  a  function  of 

PSIC 

channel  width  in  y-direction  (ny) . 

Amplitudes  of  the  cosine  function 

LAMC 

Phase  differences  for  the  cosine 

PSIS 

functions  (\pc)  in  whole  degrees. 

Amplitudes  of  the  sine  functions 

LAMS 

(^s)  in  whole  degrees. 

Phase  differences  for  the  sine 

WX 

functions (A  )  in  whole  degrees. 

o 

Adjustment  of  the  wave  near  the 

rigid  boundaries : 

Boundary,  WX  =  0; 

+1  row  from  the  boundary,  WX  =  0.33; 

+2  row  from  the  boundary,  WX  =  0.64; 
+3  and  more ,  WX  =  1 . 

86 


Fields  are  generated  according  to  the  formula: 

N 

^  -  Uy  +  H(*c)v  sin{(nx)  v(x)-j|^-(Ac)v}cos(ny)v 

+  (*s)v  sin'{(nx)v(x)-  TIs-(As)v}sin(ny)v] 


Example  Prediction  Area 


Nx  =  2  2  waves  in  x-direction  over  the  area  (dashed  line) 

Nx  =  N  N  waves  in  x-direction  over  the  area 

Ny  =  1  1  half  wave  in  y-direction  (solid  line) 

Ny  =  2  2  half  waves  in  y-direction  (dashed  line) 

Ny  =  N  N  half  waves  in  y-direction 

The  zonal  wind  speed  is  specified  with  an  arbitrary 
resolution  across  the  channel.  All  parameters  are  given  by  a 
DATA  s  t  at emen t . 
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Subroutine  STEPEXT (FI ,F2 ,F3,F4 , F5 , F6 ,F7,F8,F9 , F10 ,MY , F ,MARK , 
M,N,I1,I2,I3,I4) 

The  subroutine  computes  the  forcing  functions 


4 

E 

i=l 


mi 


4 

and  E  F 
i=l 


and  stores  the  result  in  the  fields  HM3 ,  H13  and  H23 

respectively  on  secondary  memories  (see  Flow  Diagram  B-3) . 
The  program  is  a  subroutine  to  STEP3P. 

Subroutine  Parameters  Description 


F1-F10 

MY 

F 

MARK 

11,12,13,14 


Working  field, 
y-parameter  field. 
Coriolis  parameter  field. 
Field  indicator  array. 
Computational  parameters . 

II  =  0  \V  *Vn  =  0 

X 


M,N 


11  =  1 

12  =  0 
12  =  1 
13  =  0 

13  =  1 

14  =  0 
14  =  1 
Field  vectors . 


W  -Vn  ?  0. 

X 

c • W  =  0. 

?-V\V  1  0. 

=  0  • 
9p 


03' 


ic 

3p 


¥■  o. 


|k*(|^xVo3)  =  0. 

|k*  (&Vw)  ^  o. 
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Flow  Diagram  b-3:  STEPEXT 
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Subroutine  VELPOT (KSI ,FORC,M,N, RESIDUE , ALFA) 


This  subroutine  computes  the  velocity  potential  from  a 
known  divergence  field 
V2X  =  D. 

In  finite-difference  form 


The  solution  is  performed  by  Liebmann  relaxation  with  an 
overrelaxation  coefficient  ALFA  equal  approximately  to  1.4,  but 

its  size  depends  on  the  area  and  mesh  width.  The  residual 

4*  6 

RESIDUE  must  also  be  given  (0.5*10  recommended  value) . 

Subroutine  Parameters  Description 

2D  array  for  the  x  field. 

2D  array  for  the  forcing  function. 
Overrelaxation  coefficient. 
Residual . 

Field  vectors . 

Remark:  A  first  guess  must  be  put  in  x  field  before  the 

execution. 


KSI 

FORC 

ALFA 

RESIDUE 

M,N 
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Subroutine  GRADPR(A,B,C,MARK,M,N) 


This  subroutine  computes  a  finite  difference  operation  of 
the  form  VA*?B. 


(VAVB)  ij  =  (VrtXBW'ilj  +  (B  i j 

+  (Aj+l-flj)(Bj+1-Bj>i  +  ^Aj-Aj-P  (Bj-Bj-lh 


Subroutine  Parameters 

A,B 

C 

MARK 

M,N 


Description 

Fields  for  operational  vector. 
Result  field. 

Field  indicator  array. 

Field  vectors . 
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APPENDIX  C 


PROGRAM  LISTINGS 


Three  programs  for  the  three-parameter  model  in 
Fortran  IV  are  presented:  PROG3P ,  STEP3P ,  and  STEPEXT. 
The  other  subroutines  may  be  obtained  from  ENVPREDRSCHFAC 
by  request. 
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o  o 


PROGRAM  PR0G3P 


PROGRAM  PR0G3P<0UTPUT,TAPE6*0UTPUT,DATA, lNPUT*DATA> 

BAROcLINIc  BALANCED  INTEGRATED  3-PaRaMETER  MODEL  INCLUDING  humidity 

and  precipitation,  boundary  values  can  be  variable#  constant  or  of 
channel  type, 

DIMENSION  F1(57.57),F2(57»57),F3<57.57).F4(57.57),F5(57.57) . 

1  F6(57,57),F7(57,57),F8(57,57).F9(57,57),F10(S7',57) , 

2  P( 57, 57), My (97. 57), MARK (57. 57) 


DIMENSION  UPS(10)*UPM(10)'UPl(lQ)»WX<l5)»NX<20)*NY(20)*PSIC<20)' 

X  PSJS(20),LAMC(2u),LAMS(20> 

DIMENSION  ZB<250),PSIB<290),FB<250).SbI250), 1T(2,25) 

DIMENSION  MAPHoUR(10),LETACC<10)',  NsmUTT(10),nELLIPT(10) 

DIMENSION  MSMUTT ( 10 ) 

DIMENSION  MELLIPT(IO) 

COMMON/FPAR/IC<8># JCOI.XP0L, YPOL.R.RE.DS, JM INI 100>. JMAX(IOO). 

X  M.N.KIND 

COMMON/FORMl/ 1  Cl <  8 ) , JC1 ( 8 ) , JM I N1 ( 10 0 ) . JMAX1 ( 1 0 0 ) , XP1 . YP1 . D1 
COMMON/DRM/MN.NDIM.NFLD 

common/kanal/fim 
COMMON  F 

COMMON/ECS/  PSIMl,PSIll,PSl2j.,PSIM2,PSU2,PSI22,HUMi,HUM2,DIVi,  ... 

2dIV2,WS,HEAT,J789,J12,j56, J3.PS, TS.PREC.STRM,STRl,STR2.ZMl,Zll.z2lPR0G3p,3f| 


PROG3P.2 

PR0G3P.3 

PROG3P . a 

PR0G3P.5 

PR0G3P.6 

PR0G3P.7 

PR0G3P.8 

PR0G3P.9 

PROG3P.10 

PR0G3p.li 

PR0G3P.12 

JUN12.2 

JUN12.3 

PR0G3P.13 

PR0G3P ■ 14 

PR0G3P .15 

PR0G3P.1A 

PR0G3P , 17 

PR0G3P.18 

PR0G3P ,19 

APR*  4 , 2 

PROG3P.2? 


3,ZM2,Z12,Z22,H13.H23,HM3JHM2,H12,H22,H11,H21. J4.VM.V1.V2  PR0G3P,3i 

COMMON/COEFF/Al* A2«A3.Sl.B2.B3.Cl.C2.C3»C4iC5.C6«C7.C8.D.DELP.EM.  PROG3P.32 
X  E1,E2.H1,H2,h3,h4,H5,H6,PMEAN,S1.S2,T1,T2,T3,T4,T5,  PR0G3P.33 

X  P0.PM.P1  PR0G3P.34 

COMMON/COEFF2/T6,T7,T6,T9,T1P,T11,T12»T13,T14,  PR0G3P.35 

X  K1.K2,K3.K4»K5.K6,K7«K8»K9.K10.K11.K12>K13.K14,K15*  PR0G3P . 3* 

X  K16,K17,K18,K19,K20,K21.K22,K23,K24,K25,K26,K27,K28,PROG3P.37 

X  K29,K30.K3l.K32.K33.K34.K35,K36,K37.K38  PR0G3P.38 

common/runpaR/delt,ntstep,alfasys,alfam.alfaZ.alfapsi.Rfssys,resm.prog3p,39 

X  RESZ,RESPSI,O,F0CEaN.FC0NT,WGT1,WGT2,WGT3,adIFF  PROQ3P.40 

REAL  MY  PR0G3P.41 

REAL  K1*K2.K3*K4*K5*K6»K7*K8«K9»K10*K11.K12*K13*Ki4»Ki5*  PROQ3P.43 

X  K16.K17,K18.K19,K20.K21»K22,K23.K24»K25.K26.K27.K28.PROG3P.43 

X  k29,k30,K31,k32,K33.k34,k35,k36,K37,K38  PR0G3P.44 

DATA  PSIM1,PSU1,PSI21,PSIM2,PSU2,PS122,HUMJ,HUM2,DIV1,  APR14.3 

*DI V2.WS.HEAT. J789, J12, J56. J3.PS,TS.PREC,STRM.STRi,STR2.ZMi,Zii.Z21APRi4,4 
*.ZM2.Z12.Z22.H13,H23,HM3.HM2.H12',H22.H11.H21, J4.VM.V1. V2  APR14.5 

*/l. 3250. 6499. 9748. 12997. 16246. 19495. 22744. 25993, 29242. 32491. 33740. APRI4, a 

.  TfiOfin  /  A  a  t?  a  a  **  A  O  d  m  ^  rs  fs  tr  j  ipa  i  aa  *  .  n  ^  ^  v 


*38989, 42238. 45487, 48736. 51985. 55234, 58483, 61732. 64981, 68230. 
*71479, 74728, 77977, 81226, 84475, 87724, 90973.94222,97471.100720. 
*103969, io72i8,iio467, 113716. 116965. 120214. 123463, 126712/ 
NaMELIST/FPaRR/IC.JC.XPOL.  YPOL.R'.RE.pS.JMJN.  JMAX.M.N'.KIND 
DATA(  ICKI).  I  *1. 8  )/l,  1,57. 57. 57. 97. 1.1/ 

DATA(JC1( I). 1*1, 8)/l, 1,1. 1,57, 57. 57. 57/ 

DATA  XPl.YPl.Dl  7*21. ,29., 95, 25/ 

DATA  FIM/  50,  / 

DATA  KIND/0/ 

DATA  I  CASE. I  DAY, I T I  ME/1.  660922,  0  0/ 

data  alfasys.alfam.alfaz. alfapsj/, es-,  1.45,1, 4. i.e/ 

DATA  RESSY9,Rp8M.RESZ,RE8PSI/,5E5,.5E5,.5«  ,5E5  / 

DATA  Q,FOCEAN,FCONT,w0Tl'.WGT2,  WQT3.ADIFF/  o.O  , ,  62,1 . .  1 , ,  ,5,  ,2, 
*  0,0/ 

DATA  STAB1.STAB2.PnIVS,PnIVM,PNIV1/. 422222, ,511111,100 , ,50 ,, 30 ./ 
NTSTEP  IS  NUMBER  OF  TIME  STEPS  IN  3  HOURS 
DATA  NTSTEP/4/ 

DATA  NBND/48/ 

DATA  IVAR/  6  / 


APR14.7 

APR14.8 

APRi4,9 

PROqSP ,45 

PR0G3P.4S 

PR0G3P.47 

PR0G3P , 48 

PR0G3P , 49 

PROG3P.50 

PR0G3P.51 

JUN1 2 , 4 

PR0G3P , 53 

PR0G3P.54 

PR0G3P.55 

PR0G3P.56 

APR14 • 10 

JUN12,g 

APR14.12 

PR0G3P.59 
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PR0G3P  (Continued) 


c 

c 

c 

c 


c 

c 


c 


o  aid  adiff  have  to  be  defined  *** 

MAPHOUR  GIVES  the  HOURS ( MULTIPLE  OF  6)  FOR  MaP-PR I NTI NG , 

LETACC  GIVES  THE  HOURSIHULTIPLE  OF  6)  AT  WHICH  THE  ACCUMULATED 
SHALL  BE  PUT  =0  AGAIN. 

MSMilTT  GIVES  THE  HOURS < MULT  I PLE  OF  6>  FOR  SMOOTHING, 
f  ELL  I PT  GIVES  The  HOURS  (MULTIPLE  Op  6>  FOR  ELLIPTIcITY  TEST. 
DATA  MAPHQIIR/0,12,24,3*,48,-1,-1.-1.-1*-1/ 

Data  (LETaCC<I>,I=1»10)Z0»12. 24,36, 48,-1, *1,-1, -1,-1/ 
DATA(NSMUTT( I). 1  =  1, 10  I/O. 3, 6. 9, 12. 15. 18,21,24,-1/ 

DATA  (MSMUTTf I)» 1  =  1*10 1/27* 30* 33*36* 39 *42*  45*48, -1*-1/ 

DATA (HELL  I  PTC  I ), I =1, 10  I/O , 3, 6, 9, 12, 15, 18, 21, 24, .1/ 

DATA  (MELLIPT( I ). 1=1*14 )/27, 30.33, 36,39*  42, 45, 48,-1, -1/ 
3=9,806 


PR0G3P , 6  0 
PR0G3P.61 
PPECIPIPR0G3P.6? 

PR0G3P , 63 

PR0G3P .64 

PROG3P.65 

JUN12.6 

JUN12.7 

PR0G3P.6A 

JUN12.8 

PROG3P , 69 

JUN12.9 

PROG3P.70 


F0*1.03E-4 

TIME  STEP  IN  SECONDS 
UEUT=i.*3.6E3/NTSTEP 
LABEL  =  lOH  PUTl 
3 T  I ,kI E  =  SECONpCFAKE) 

WRITE(6*8000>  BTIME 
eUOO  FORMAT ( lHO ,  *BTIME=*,  F1H.4) 

CALL  RUT1 

OTIME  =  SECONd(FAKE)  -  B  T I  HE 
WRITE(6*H005>  LABEL*  DTImE 
8n0s  FORMAT ( lHO ,  -TIME  TO  EXECUTE*.  AlO.  F10.4) 
LABEL  =  10H  COEFF3P 
3TIME  =  SECOMU(FAKE) 

WRI TE ( 6, 8000 )  BTIME 

CALL  C0EFF3P(STAB1.STAR2.PNIVS,PNIVM.PNIV1) 

DTIME  =  SECOND(FAKE)  -  BTIME 

-IR  I  TE  (  6, 8°  05  )  LABEL,  DTIME 

MN  *  M*N 

NDIM  =  20000 

LABEL  =  10H  NaRKF 


PR0G3P . 7 1 

APR14.14 

JUN12.10 

PR0G3R.73 

PR0G3P.74 

PR0G3P . 75 

PR0G3P.76 

PROG3P.77 

PR0G3P.7P 

PROG3P . 79 

PROG3P.80 

PR0G3P.81. 

PR0G3P.8? 

PR0G3P ,83 

PR0G3P . 84 

PR0G3P.85 

PR0G3P , 86 

PR0G3P . 87 
PR0G3P . 88 
PR0G3P.89 


BTIME  =  SECOND(FAKE) 
WRITE(6,8000)  BTIME 
CALL  MARKF(MARK,M,N) 

DTIME  =  SECOND(FAKE)  -  BTIME 
WR I TE ( 6, 80 05 )  LABEL,  DTIME 
LABEL  =  10H  MYFF 
BTIME  =  SECOND(FaKE) 
WRITE(6,8°00>  BTIME 
CALL  MYFF(MY.M.N) 

DTIME  =  SECONd(FaKE)  -.BTIME 
WR I TE ( 6 » 80 05 >  LABEL*  DTIME 
DO  9  1=1, MN 
9  FI (  I  >  =  0  ,  0 

LABEL  =  10H  RANWT 
BTIME  =  SECOND(FAKE) 

WRITE ( 6, 8000 )  BTIME 
CALL  RANWT(VM.Fl) 

CALL  RANWT ( VI, FI ) 


PR0G3P . 90 
PR0G3P , 91 
PR0G3P,9? 
PR0G3p,93 
PR0G3P .94 
PR0G3P , 95 
PR0G3P.96 
PR0G3P.97 
APRl4,l6 

PR0G3P,9« 

PROG3P.100 

PR0G3p.l*)l 

PR0G3p,m2 

pROG3P,103 

PR0G3P.1Q4 

PROG3P.105 

PROG3P  .1(|6 

PR0Q3P , 107 


CALL  RANWT ( V2. FI > 

DTIME  =  SECOND(FAKE)  -  8T I ME 
WRITE (6, 8005)  LABEL,  DTImE 


PR0G3P . 108 
PROG3P.109 
PR°G3PtliO 


NTIME  =  0 

IF(KINn.NE.O)  GO  TO  45 
C  INITIAL  FIELDS,  NO  CHANNEL. 

C  THE  subroutine  INITF  HAS  TO  BE  WRITTEN  *** 

C  INITIAL  Z-FIELDS  at  the  LEVELS  PS. PM  and  P1  TO  F*  *  F2  and  F3.  the  HUMIDPROG3P. 
C  FIELD  Q  TO  F 4.  STANDARD  SURFACE  PRESSURE  PS  TO  F5  AND 


PR0G3P ,  111 
PR0G3P.U2 
PR0G3P ,113 
PR0G3P ,114 
"  115 

SEA  LEVEL  TEMPERPROG3P . 1$ 6 
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PR0G3P  (Continued) 


C  TS  TO  F6,  DEFINITION  Q*0, 0334892*M ( 850 > *0 , 0363ft89*M< 700 > *0 . 029095l*M C PR0G3P . U7 
C  0,0145476*M<300>  WHERE  M(r>  ARE  MIXINQ  R*TI°S  jN  MTS»UnITS.  PS*101 . 325PR0G3P , 110 


C  AT  SEA  SURFACE. 

LABEL  ■  10H  INITF 
BTIME  ■  SECOND(FAKE) 

WR I TE ( 6, 8® ® ®  >  BTIME 

CALL  INITF<Fl.FZ.F3,F4tF8,F6.F7(l.l),F7(l,5) ,F7(l,9),F7(l,l3), 

*F7(l,l7),F7(l,21),F7(1.25),F8fM,N) 

SCALE*20. 

DTIME  a  SECOND(FAKE)  •  BTIME 
WRITE<6.80°5>  LABEL.  DTIME 
LABEL  a  10H  MAP  6X 
BTIME  a  SECOND(FaKE) 

WR I TE ( 6, 8®  ® ®  >  BTIME 
C  SOLUTION  of  streamfunctions 
LABEL  a  10H  STREAMF  1 
BTIME  =  SECOND(FAKE) 

WRITE(8.8000)  BTIME 
LABEL  »  10H  STREAMF 
BTIME  a  SECOMD(FAKE) 

WRITE(6»8000>  BTIME 

Is  call  STREaMF<F1»F7»F1®.MaRK,ZB.PSIB.FB.SB.GAMMA,0, fTl.M.N) 


PR0G3P.119 

PROG3P.120 

PR0G3P.121 

PROG3p.l?2 

PR0G3P.123 

APR14.17 

PR0G3P.125 

PR0G3P , 1?6 

PROG3P.1P7 

PR0G3P . 128 

PR0G3P . 129 

PR0G3P.13® 

PR0G3P.139 

PROG3P.140 

PROG3P .141 

PR0G3P ,142 

PR0G3p,l43 

pR°G3P,144 

PR0G3P .145 

APR14.18 


DTIME  =  SECOND(FAKE)  -  BTIME 
4RITE(6.8005)  LABEL,  DTImE 
LABEL  *  10H  STREAMF  2 
BTIME  a  SECOnD(FAKE) 

WR I TE ( 6 •  80  0  0  )  BTIME 

CALL  STREAMF(F2,F8,F10,MARK,7B,PSIP,FB.SB,GAMMA,0, IT2.M.N) 

DTIME  =  SECOND(FAKE)  -  BT I  ME 

WR I TE ( 6, 80 05 )  LABEL,  DTIME 

LA8EL  =  10H  STREAMF  3 

BTIME  =  SECONDCFAKE) 

WRITE(6,8000)  BTIME 

CALL  STREAMF C F3,  F9,  FI  0.  MaRK,  7B.PS I B.FB.  SB.  GaMM A,  0,  IT3.M.N) 

DTIME  =  SECONDIFAKE)  -  BTIME 

WR I TE (  6, 80 05 )  LABEL,  DTImE 

DO  20  1*1, MN 

F5( I )=F5< I ) *0 , 1 

Z  =  F2(  I  ) 

F2< I >*.5*(Z«F1«  I  )  ) 

F3( I >=.5*<F3< I >-7> 

Fl<I)*Z 

20  CONTINUE 

C  ADAPTION  OF  HUMIDITY  FIELD 
DO  21  1*1, MN 

TZ  =  G*(H3*F2< I )*H4*F3( n )/F0 

TPS  I  =  , 5*(H3*(F0(I)„F7(!))+H4*(F9(I).F8(I>)) 

CALL  SATUR(TZ.OZ) 

CALL  SaTUR(TPSI.OPSI) 

C  PHaNgE  IN  F 4  SINCE  F4  IS  RELATIVE  HUMIDITY 
F4< I )*F4< I ) *OpS I 
C  F4( I )*F4< I J*DRSJ/QZ 

Z  =  F4( I >-.8*QPSI 
IF(Z.GT.O.O)  f4< I )=.8*0PSI 

21  CONTINUE 

C  STORE  Z-FIELDS  FUR  BOUNDARY  MIXING 
LABEL  =  10H  RANWT  3X 
BTIME  =  SECOND(FaKE) 

WRITE(5,8t)00)  BTIME 
CALL  RANWT(ZM2.F1) 

CALL  RaNhT ( Z12. F2 ) 


PR0G3P .147 

PROG3P.148 

PR0G3P.149 

PR°G3P.150 

PR0G3P  ,151 

APR14.19 

PR0G3p,153 

PR0G3P.154 

PR0G3P . 155 

PR0G3P.156 

PR0G3P , 157 

APR1 4,20 

PR0G3P , 159 

PROG3P.160 

PR0G3P.1A1 

PROG3p,ls2 

PR0G3P.163 

PR0G3P.1A4 

PR0G3P.165 

PR0G3P.166 

PR0G3P.187 

PR0G3P .188 

PROG3P . 1*9 

PROG3P.170 

PR0G3P. 171 

PR0G3P.172 

PR0G3P,173 

PR0G3P,‘l74 
PR0G3P.175 
PR0G3P . I76 
PR0G3P , I77 
PR0G3P ,178 
PROG3P , 179 
PROG3P.I8Q 
PROG3P ,181 
PR0G3P , la2 
PR0G3P.183 
PR0G3p , lfl4 
PROgjP ,185 
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PR0G3P  (Continued) 


CALL  RANWT ( Z22*F3 ) 

PR0G3P.186 

DTIME  s  SECOND(FAKE)  •  BTIME 

PR0Q3P.187 

WR I TE ( 6 , sOOg )  LABEL,  DTIME 

PR0G3P.188 

test  for  ellipticity  of  psi-fields 

PR0G3P.l89 

LABEL  ■  10H  ELLIPT  3X 

PROG3P.190 

BTIME  »  SECOND(FAKE) 

PR0G3P ■ 191 

WRITE(6,8000)  BTIME 

PROG3P , I92 

CALL  ELLIPTCF7,MY,M,N> 

APR14.21 

CALL  ELL  I PT ( F8. MY, M , N ) 

APR14.2Z 

CALL  ELLIPT(F9,MY,M,N) 

APR14.23 

DTIME  »  SECOND(FAKE)  -  BTIME 

PR0G3P.196 

WR  I  TE  ( 6, 80  O5  >  LABEL,  DTIME 

PROG3P . 197 

computation  of  psim.psu  and  psi? 

PROG3p.190 

DO  25  1*1, MN 

PR0G3P.199 

Z=FS< I ) 

PROG3P.2QO 

F8 1 1 >«,5*(Z«F7<I ) ) 

PR0G3P .201 

F9( I )  =  .5*(F9( I > -Z  > 

PROG3P , 2o  2 

F7< I )*Z 

PR0G3P , 20  3 

25  CONTINUE 

PROG3P.204 

IF( IVAR.EQ.O*  GO  TO  31 

PROG3P .  2fl5 

STORE  FIELDS  FOR  BOUNDARY  “IXlNG 

PR0G3P ■ 206 

DO  30  1=1, MN 

PR0G3P .207 

Z  =  F1<  I  ) 

PR0G3p.2fl8 

F 1 C 1 1 =F7 (  I  ) 

PROG3P.209 

F7< I >=F7< I ).G*Z/F0 

PR0G3P ,  2i 0 

Z=F2< l\ 

PR0G3P.2H 

F2( I >=F8<  I  ) 

PROG3P , 212 

F8< I >=F8< I )-G*Z/F0 

PR0G3P . 213 

Z  =  F3(  1  ) 

PR0G3P.214 

F3(  I  )=F9< I ) 

PR0G3P.215 

F9< I )=F9( I >-G*Z/FQ 

PR0G3P . 21 6 

30  CONTINUE 

PR0G3P.217 

GO  TO  36 

PR0G3P . 21 8 

31  DO  35  1=1, MN 

PROG3P . 21 9 

F1C I >=F7(  I ) 

PROG3P,2?0 

F2(  I  )=F8( I ) 

PR°G3P,221 

F3(  I  >=F9( I > 

PR0G3P , 272 

35  CONTINUE 

PR0G3P.273 

LABEL  =  10H  RANWT  3X 

PR0G3P .  2?4 

3TIME  =  SECOND(FAKE) 

PR0G3P , 2?5 

WR I T£(6, 8000 )  BTIME 

PROG3P.226 

36  CALL  RANWT ( STRM, F7 > 

PR0G3P.227 

CALL  RANWT  t  STR1, F8 ) 

PROG3P.2?8 

CALL  RaNwT(STR2,F9) 

PR0G3P .279 

DTIME  =  SECOND I FaKF  )  -  8 T I ME 

PR0G3P . 230 

WRITE16.8005)  LABEL,  DTIME 

PROG3P.231 

PRINT  NUMBER  OF  ITERATIONS  in  THF  SOLUTION  OF 

psi 

PROG3P. 232 

41  WRITE(6,202)  IT1.IT2.IT3 

PR0G3P . 233 

p0?  FORHaT(///1Xi33HNUM0ER  OF  ITERATIONS  IN  STREAMF  A3d4)) 

PR0G3P.2,34 

GO  TO  53 

pR°G3p.235 

generation  of  initial  fields  for  the  CHANNEL 

VfRSION,  PSI-FIELDS  aPE 

PROG3P . 236 

generated  in  the  subroutine  bench  from  Parameters  given  in  daTa-staTemprog3P.237 

HUMIDITY-FIELD  CORRESPONDS  to  50  PERCENT  RELATIVE  HUMIDITY.  PS  IS  loo  PROG3P.238 
AND  TS  273  DEGREES  EVERYWHERE.  PR0G3P.239 

45  CALL  GENCH < F7. M, N, PS1 PS, UPS. NU.NWaVE. NX. NY. PS JC.LaMC.PS IS, LAMS.WXIPR0G3P. 240 

WRITE(*.FPARR)  PR0G3P.241 

call  GEMcH ( F8, M, N, Ps I  pH,, jPm, Mu, nKA VE, Nx, ^Y, PS  I C.LAMC.PS IS, LAMS, WX)PR°G3P. 24  2 
call  GENCH(F9»M»N*PSIPl»UPi»NU»NWAVE*NX»NYiPSIC»UMC*PSIS»LAMS.WX)PR0G3P.243 

COMPUTATION  OF  PSIM.PSI1.PSI2.HUMIDITY.PS  AND  TS  PROG3P.244 

46  OO  5q  1=1. MN  PR0G3P.245 
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PR0G3P  (Continued) 


FlC I )=F8<  i  > 

F2< I >  =  ,5*(F8< I >*F7( I >  > 

F3<l)«,5*(F9(I)"F8<n) 

F5(J)»  loo. 

F6< I >  =  273. 

TPSIs  H3*F2( I )*H4*F3(  I  ) 

CAUL  satur‘tpsi,qpsi > 

F4( I )=.5*QPSI 
50  CONTINUE 

c  STORE  PSI-FIELDS  In  STRM.STRi  aNd  STR2 
CALL  RANWT(STRM*F1> 

CALL  RANNTCSTRl.F2) 

Call  RANWTCSTR2,F3> 

C  SMOOTHING  OF  INITIAL  fields 

53  I r < NSMUTT < i > . NE , 0 )  GO  To  55 
LABEL  =  10H  ASMUT  6X 
3TIME  =>  SECOND!  F  AKE ) 

WRITEC6.8000)  BTIME 
CALL  A9rtUT(Fl.FlO.M,N»  .5) 

CALL  ASMiJTfFl.FlO.M,  N,„(g) 

CALL  ASMUT(F2»FX0*M*N»  .5) 

CALL  ASHUT(F2,FlO#M,N,-,5> 

CALL  ASMUTCF3,F10,M,N,  .5) 

CALL  ASMUT(F3,F10,M,N,.i5) 

DTIME  =  SECOND(FaKE)  -  BTIme 
RRITE(6.8005)  LABEL.  DTImE 

C  LOADING  OF  INITIAL  FIELDS.  ZcRO  TO  ACCUMULATED  PRECIPITATION 

55  DO  56  1=1, MN 
F10(I)»0.0 

56  CONTINUE 

LABEL  *  1 0 H  RANWT  10X 
BTIME  =  SECONDCFaKE) 

WRITE(6,8000)  BTIME 
CALL  RANWTCPSIMl.Fl  ) 

CALL  RANWTCPSIU.F2  ) 

CALL  RANWTCPSI21.F3  ) 

CALL  RANWTCHUM1,  FA  ) 

CALL  RANWTCPS  ,F5  ) 

CALL  RANWTCTS  . F6  ) 

CAUL  RANWTCPREC  , F10 ) 

CALL  RANWTCWS  ,F10) 
call  RANWT ( D I VI  « Fl 0  ) 

CALL  RANWTtDIV?  .FlO) 

DTIME  =  SECONDCFaKE)  -  BTIME 
WRITEC6.8005)  LABEL,  DTIME 
IFCMaPHOURCD.NE.O)  GO  TO  60 
C  MAPPRINTING  of  FIRST  TIMESTEP 
IFCKIND.NE.O)  go  TO  58 
LABEL  =  10H  MAP3P 
3TIME  =  SECOMD(FAKE) 

WRITEC6.8000)  BTIME 

CALL  MAP3P(1, 1,1,1,0»2,0.2, 0,0,0. FI, F2,F3,F4,F5,F6,F7,F8,F9, Fin, 
X  MARK.M.N, I  DAY, IT IME, NjIME.PN I  VS, PN iv M, PNIvi.78.PS IB,  FB) 

DTIME  =  SECONDCFAKE)  -  BTIME 
WRITEC6.8005)  LABEL,  DTIME 
GO  TO  60 

58  CALL  MAP3P(0,0,0,0.0,0,0.0.1,1,1,Fi,F2,F3,F4.F5,F6,F7.F8,F9,F10, 
X  MARK.M.N. I  day, ITIME,NTIME,PNIVS,PNIVM,PnIV1.2B«ZB,ZB) 

LABEL  =  10H  FORECASTS 
BTIME  =  SECONDCFaKE) 


PROgSP , 246 
PR0G3P .247 
PR0G3P , 248 
PR0G3P.249 

PROG3P.250 
PRUG3P.251 
PR0G3P ,252 
PR0G3P .253 
PR0G3P.25« 
pR°G3P.255 
PR0G3P , 256 
PR0G3P.257 
PR0G3P,258 
PR0G3P • 259 
PR0G3P • 260 
PR0G3P.261 
PR0G3P.262 
PROG3P.263 
PR0G3p ,264 
PR0G3P.265 
PR0G3P .266 
PR0G3P.267 
PR0G3P , 268 
pR°G3P.269 
PR0G3P . 270 
PR0G3P . 271 
PROG3P.272 


PROG3P . 273 
PROG3P , 274 
PR0G3P.275 
PR0G3P , 276 
PR0G3P , 277 
PROG3P , 278 
PR0G3P . 279 
PROG3P.280 
PROG3P.281 
PR0G3P  ,282 
PR0G3P ,283 
PR0G3P ,284 

PR0G3P .  2a5 

PR0G3P.286 
PR0G3P.287 
PROG3P.208 
PR0G3P.289 
PR0G3P . 29O 
PR0G3P .29I 
PROG3P.292 
PR0G3P . 293 
PR0G3P.294 
PR0G3P.295 
PR0G3P . 296 
APR14 ,24 
PR0G3P . 298 


PR0G3P .29? 

PROG3P.300 

PROG3p,3ol 

APR14.25 

PR0G3P .303 

PROG3P.304 

PROG3P , 305 
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PR0G3P  (Continued) 


WR I TB ( 6, 8000 )  BTIME 
C  THREE  HOURS  FORECAST  STARTS  HERE 

60  DO  61  1*1, MN 

F1U>  =  0,0 

61  CONTINUE  1 

CALL.  RANWT<DIV1,F1> 

CALL  RANWT ( D! V2«  FI ) 

CALL.RANRD(HUM1,f2) 

CALL  RANwT(HUM2,r2) 

CALL  RANRD ( PS  I  Ml*  F2 ) 

CALL  RANWT ( PS  I  M2  , F2) 

CALL  RANRD(PSJ11,F2) 

CALL  RANWTI PS  1 12, F2 ) 

CALL  RANRDlPsI21,F2) 

CALL  RANWT ( PS  1 22, F2  > 

IFUVAR.EQ.0)  QO  TO  70 
IF(KiNd.Ne.o)  QO  TO  70 
C  INPUT  OF  BOUNDARY  FIELD  EACH  SIX  HOUR 
C  THE  SUBROUTINE  ZINPUT  HAS  TO  8E  WRITTEN  *** 

CALL  RANRD(ZM2,F2) 

CALL  RANWT(ZM1,F2) 

CALL  RANRD(Z12.F2) 

CALL  RANWT(Z11,F2) 

CALL  RANRD ( Z22, F2 ) 

CALL  RANWTIZ21.F2) 

C  CALL  ZINPUTO  ZPS,ZPM,7P1  TO  FIELDS  F1.F2  AND  F3, 

DO  65  1*1, MN 
Z  =  F2(  I  ) 

F2< I  )i ,5*(Z-F1( I ) ) 

F3( I >=.5*(F3( I )-Z) 

Fl< I )=Z 
65  CONTINUE 

CALL  RANWT<ZH2.Fl ) 

CALL  RANWTIZ12.F2) 

CALL  RANWT ( Z22*  F3 ! 

C  GENERAL  THESTEP  THREE  HOURS  AHEAD 

70  CALL  STEP3P(Fi,F2,F3,F4,F5.F6,F7.F8,F9.F10.My.MARK, IT, IvAR.H.N, 
NCOUNT  *NT I  ME 

NTIME*NTIME+1 

u=ntstep*i 

C  PRINT  NUMBER  of  ITERATIONS 

WRITE(6,200)  NCOUNT, NTIMF 
DO  71  1*1,11 

HR  I TE ( 6, 201 )  IT(1,I),IT(2,I) 

71  CONTINUE 

C  SMOOTHING  AND  ellipticity  test 
11  =  0 
12  =  0 

DO  80  Isl,10 

IF (ABS( FLOAT (NSMUTTt I > ) -NTJME ) ,LT . O ,1 )  11=1 
IF(ABS(FLOat(MSMUTT( I n.NTIME) ,LT.0,1>  11=1 
IF< ABS(FLOAT(MELLIPT(I))«NTIME).LT.0.1>  12=1 

IF<ABS<FLOATCNELLIPT<1 ) > -NT  I  ME) .LT.0,1)  I  2*1 
80  CONTINUE 

IFUl.EQ.O.AND.  l2,EQ.Q)  GO  TO  95 
CALL  RANPD ( PS  I  Ml, FI ) 

CALL  RANRD(PSIU.F2) 
call  RANRD  < PS  1 21 , F3 ) 

DO  85  lal.MN 
Z  =  F1 (  I  )  • 


PR0Q3P.3&6 

PROG3P.307 


PROGJP , 3q8 
PROG’P ,  *0  9 
PR°G3P,3i0 
PR0G3P.311 
PR0G3P .  3l  2 
PROQSp  3«3 
PR0Q3Pi3l4 
PR0G3P , 3{5 
PR0G3P , 316 
PR0G3p,3i7 
PR0G3P ,318 
PR0G3p,3l9 
PROG3P.320 
progSp,3?i 
PR°G3P,322 
PR0G3P . 3?3 


PR0G3P.324 
PR0G3P.325 
PR0Q3P ,356 
PR0G3P .  3?7 
PROG3P.328 
PR0G3P , 329 


PR0G3P , 330 

PR0G3P.331 

PR0G3P.332 

PROG3P .333 

PR0G3P.334 

PR0G3P.335 

PROG^P , 336 

PR0G3P,337 

PR0G3P.338 

PR0G3P.339 

PR0G3P , 340 

PR0G3P.341 

APR14.26 

PR0G3P.343 

JUN12.il 

PR0G3P.345 

PR0G3P , 346 


PR0G3P.347 

PR0G3P.348 

PR0G3P.349 


PROG3P .350 
PR0G3P ■ 351 
PR0G3P , 352 
PR0G3P ,353 
PR0G3P.354 
JUNl2,l2 
JUN12,l3 
JUN12.14 
JUN12.15 
PR0G3P , 357 
PR0G3P.358 
PR0G3P , 359 
PROG3P,360 
PR0G3P , 3*1 
pR0G3p ,362 
PR0G3P.363 
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PR0G3P  (Continued) 


F*<  I  )  =  Z«2*F2( i ) 
f3(D  «  z*2.*r3m 

F2(  I  )  =  Z 

85  CONTINUE 

I F  (  U.EQ.O)  GO  TO  86 
CAUL  A3MUTJF1.F4.M.N,  ,5) 

CALL  ASMUT(F1.f4»M.N.".5) 

CALL  ASMUT(F2,F4,M,N,  ,5, 

CALL  ASMUT<F2'F4»M»N*->,5) 

CALL  ASM'JT  ( F3,  F4 ,  M,  N,  ,5) 

CALL  ASMUTtF3,F4.M.N.-.5>' 

86  IF< I2.EQ. 0)  GO  TO  90 
IF(KiNd,NEi0>  qO  TO  90 
CALL  ELLIPT(Fl*MY>M.N) 

CALL  ELLIPT(F2,MY,M,N> 

CALL  ELLIPT(F3,MY,M,N) 

90  DO  93  1*1, MN 
Z  =  F2(  I  ) 

F2< I )s,5*(Z«Fl(  I  )  ) 

F3<  I)*.5*(f3(  n-Z> 

Fl( I )*Z 
93  CONTINUE 

CALL  RANWT ( PS  I  Ml, Fl ) 

CALL  RaNWT ( PS  1 11 . F= ) 

CALL  RANWT < PS I21»F3) 

C  MAPPRINTING 

95  DO  96  1=1,10 

IFCMaPHOURI I ) .EQ.NTIME'  GO  TO  97 

96  CONTINUE 
30  TO  100 

97  IF(KIND.NE.O)  GO  TO  98 

CALL  MAP3PI1, 1,1, 0,1, 1,1,1, 0.0. 0,Fl,F2«F3,F4,F5,F6,F7,F8,F9,Fl0 
X  MARK.M.N, I  day, ITIMF.NTINE.PNIVS.PNIVM.PNIVI.ZB.FB.SB) 

GO  TO  100 

98  CALL  MAP3P(0,0,0,0.1,0.0,0.1.1,0,Fl,F2,F3,F4,F5,f  6,F7,F8,F9,F10. 

X  MARK* M »N» IDA Y» ITIMEf NTlMt*PNIVS*PN!VM»PNIVl»7B»Z8*Z8) 

C  7ER0  TO  ACCUMULATED  PRECIPITATION 
100  DO  105  1=1.10 

IF(LETaCC( I ) .EQ.NTIME)  GO  TO  106 

105  CONTINUE 
GO  TO  115 

106  JO  110  1=1. MN 
Fl (  I  )  =  0 , 0 

110  CONTINUE 

CALL  RANWT(PREC.Fl) 

200  FORMAT <29HlNUM9ER  OF  ITERATIONS  HETWEEN. 1 3 . 4H  AMD.I3.6H  HOURS) 

?0l  FORMAT(4x,2( I  4  > ) 

115  IF(N0No,LE,NTIME)  STOP 

JU  To  60 
END 


PR0G3P , 3ft4 
PR0G3P,345 
PR0G3p,366 
PR0G3P.367 
PR0G3P.368 
PROG3P , 3$9 
PROG3p,370 
PROG3P ,371 
PR0G3P , 372 
PR0G3P ,373 
PROG3P , 374 
PR0G3p , 375 
PR0G3P.376 
APR14.27 
APR14.28 
APR14  29 
pROG3* ,  38  0 
PR0G3P .381 
PR0G3p.3fl2 
PR0G3p.3fl3 
PR0G3P , 3p4 
PR0G3P , 385 
PR0G3P.386 
PR0G3p,3fl7 
PR0G3P.388 

PR0G3P.3a9 
PR0G3P , 390 
PR0G3P.39l 
PR0G3p , 392 
PR0G3P.393 
PR0G3P.394 
APR14 ,30 
PROG3P ,  396 
PR0G3p,397 
APR14.31 
PR0G3P.399 
PR0G3P , 400 
PROG3P.401 
PR0G3P , 4  n2 
PROG3P.403 
PR0G3P , 4  0  4 
PR0G3P . 4fl5 
PR0G3p.4n6 
PR0G3P , 4  0  7 
PR0G3P , 408 • 
PROG3P  .  4n9 
PROG3P.4l0 
PROG3P ,  4H 
PROG3P , 4{2 
PR0G3P.413 
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SUBROUTINE  STEP3P 


SUBROUTINE  STEP3P<F1,F2,F3,F4,F5,F6.F7.F8,F9,F10,My,MARK, IT, 

X  IVAR.M.N) 

C  SIX  HOURS  FORECAST  BY  3p. MODEL 

DIMENSION  Fi<i>,F2<1>«F3<1),Fa<i>,F5<1>.F6<1>»F7(1>.F8(1>.F9(1>, 
X  F1D(1),MY(1),F(1).MARKU),  IT(2.1) 

C0MM0N/FPAR/IC(8), JC(8),xPOL, YPOL.R.RE.DS, JMIN(lOO). JMAX(IOO) , 

X  MX. NX, KIND 

COMMON /ECS/  PS  I  Ml , PS  1 11 » PS  1 21 , PS  I  M2, PS  1 12 , PS  1 22 , HUMj, HUM2 » D I Vj , 


APR14.46 
STEP3P.3 
STEP3P.4 
STEP3P.5 
STEP3P.6 
STEP3P.7 
STEP3P , 8 
STEP3P.1P 


201  V2.WS.HEAT. J789* J12, J56, J3.PS, TS. PREC . STRM . STR1 , STR2 . ZM1, Zll . Z21STEP3P . 19 

STEP3P , 2n 


3.ZM2,zl2,z22.Hl3,H23,HM3;HM2,H12,H22,Hn,H21. J4, vM.vl.V2 
COMMON/COEFF/A1* A2»A3>Bl,B2.B3»Cl»C2.C3«C4,C5.C6.C7»C8.D.DEl.P.EM.  STEP3P|2i 
X  E1.E2,H1,H2,h3,h4,H5.H6,PMEAN,s1,S2,T1,T2,T3,t4.T5.  STEP3p..22 

X  PO , PM, PI 

COMMON/COEFF2/T6.T7.T8, T9, Tip , Til , T12, T13, T14 , 

X  Kl,K2»K3.K4.K5,K6.K7,K8.K9,KlO.Kll,K12«K13.K14,Kl5> 

X  K16,Kl7,Kl8,Kl9,K20,K21,K22,K23.K24.K25,K26.K27,K28,STEP3P.2A 

X  K29,K30,K3l,K32,K33,K34,K35,K36,K37,K38  STEP3P.27 

COMMON/RiJNPAR/DELT.NTSTEP. ALFASYS,AUFAM,ALFAZ,ALFAPSI,RESSYS.RESM,STEP3P.2R 


STEP3p,23 

STEP3P.24 

STEP3P.25 


COMMON  F 
REAL  MY.KEFF 
REAL 


RESZ»RESPSl,Q,FOCEAN,FCONT,WQTl.WGT2.WGT3. ADIFF 


STEp3p , 29 
APR14.X7 

.  ,  „  STEP3P . 38 
K1,K2,k3,k4,k5,K6,K7,K8.K9,K10.K11,K12.K13,K14,K15.  STEP3P.3i 
K16,Kl7,Kl8,Kl9,K20,K21,K22,K23,K24,K25,K26,K27,K28,STEP3P  3? 


X  K29.K30.K3J,»K32.K33.K34.K35.K36,K37,K38 

DATA  RfiAS,EE.HL.CP.To.PO.DELl,DEL2,TOL/287  ,, .622. 2.5E6. 1004 , , 
X  273,, ,611,0,0.0.0.0,0/ 

C  DEL1.0EL2,  A JD  TUL  HAVE  TO  BE  DEFINED  ******* 

DEL1=0 . 0001 

DEL2=0.001 

TOL=0,0 

c  constants  for  computation  of  latent  heat 

CCl  =  l./TO 


CC2 

CC3 

CC4 

cc5 


EE*HL/RQAS 

EE*HL/CP 

CC2*CC3 

1EL1+DEL2 


STEP3P.33 
STEP3P.34 
STEP3?,3? 
STEP3P.36 
STEP3P.37 
STEP3P.3B 
STEP3P.39 
STEP3P.40 
STEP3P.41 
STEP3P.4? 
STEP3p , 43 
STEP3P , 44 
STEP3P.48 


c 

STEP3P.4* 

lrtKINO.EO. 0)  ou  TO  5 

STFPXP .47 

H  0  T 1  =  1 , 0 

STEP3P, 4" 

HGT2= , 67 

JUN12.16 

H3T3= . 33 

JUN12 ,17 

5 

mn=m*n 

STEP3P . 51 

'iO  =  NTSTEP  +  1 

step3p.5? 

M1=M-1 

STEP3P , 53 

DO  170  KTsl,  ID 

STEP3P . 54 

EPS  =  2, 

STEP3P .55 

I F  <  KT , LT . 3 )  EPS  =  ,5*KT 

STEP3P  ,5ft 

C 

STEP3P.57 

C  *  *  *  ★ 

**************  ********** jiCOB I  AN  COMPUTATIONS***** 

***»*♦**.«****. **STEP3P. 58 

C  all 

JACOrIaNS  are  COMPUTED  AND  STORED 

STEP3P ,  5<}' 

CALL  RANRDIPSiMi.FD 

STEP3P , 6n 

CALL  RANRD(PSI11,F2) 

STEP3P . 6l 

CALL  RANRDIPSI21.F3) 

STEP3P , 6? 

c 

STEP3p . 63 

CALL  ABSV0R(Fi,F4.MY.MARK.M,N) 

APR14 ■ 48 

CALL  RELV0r({F2,F5,MY,MARK,M,N) 

APR14.49 

CALL  RELV0rt(F3.F6.HY,MARK,M,N) 

APR14 ,50 

keff«i,  o' 

STEP3P , 67 

CALL  RaNRD(PS.MY) 

STEP3P.68 

CALL  0 »&tt<F5*MY.M.N.KEFF> 

STEP3P . 69 

r 

STEP3P  ,  70 
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STEP 3 P  (Continued) 


DO  10  I  =1 , MN 

T7(J)  ii  c3*T4(  I).C2*F'5J  I)*c4*F6(  I  )  *  C7*F<I) 
F0<!>  ■  »C2*F4(n*C5*F5<I) 

F9d>  «  04*F4(I)  ♦C6*F6(  I  )*2*C7*F<  I  > 

*0  CONTINUE 

CALL  JACOB(Fl,F7,Fl0,M,N> 

CALL  JACOB<F2,F8,F7,M,N) 

CALL  JAC0B(F3.F9.F8,M,N) 

CALL  0R0g<F7.My,M,N,kEFF> 


DO  20  1=1, MN 

F10< I)  =  F10 ( I >+F8{ [>+F7( I ) 

F9  (  I  )  a  F4  (  n.2*F5C  ) 

20  CONTINUE 

CALL  RANWT ( J789, F10 ) 

CALL  JACOB (F2»F9*F7»N«N) 

CALL  OROG(F7,MY,M,N,KEFF) 

CALL  JACO0(F1,F5.F0,M,N, 

00  30  1  =  1,, 'IN 

FI 0  C  I  )  a  F7( I ) *FB ( I ) 

F9< I )  =  Fa ( I )  +2*F6 ( I > 

30  CONTINUE 

CALL  RANWT< J12.F10) 

CALL  JACOB ( F3 , F9, F7  ,M,N) 

CALL  JACOB ( Fl »  F6»  F8  *H»N) 

CALL  JAC0B(Fl,F3,F9  ,  M ,  N ) 

CALL  JACUB(F1,f2,F10,M,N) 

CALL  OflOG(Fl0,MY,M,N,KEFF) 

CALL  RANWT( J3.F10) 

DO  40  1=1, HN 
F10 ( I ) =F7 ( I ) *F8 ( I  ) 

CALL  RAnWT( J56.F10) 

CALL  RANRD(PS«F10  ) 

ppM=2,/(p0.pH) 

DO  44  1=1, MN 

F7( I )=F1( I )-f2( I )*RPM*(F10 ( I )»PN> 

F8( I )=PPM*F5( I >  * ( FlO { I ) -PM) -F4 ( I)*F( I) 

44  F 4 (  I  )=F10<  I  ) 

CALL  JAC0B(F7.F4»F5.M, V) 

CALL  MYFF ( MY , M«  N ) 

r***************************^0^(:R  BO{JNDaRy  CONDIT 
1  NFL JENCE  FROM  TOPOGRAPHY  AND  FRICTIO  OVERR  LAND 
OCEAN  SURFACE  IS  ASSUMED  WHERE  STANDARD  PRESSURE 
DO  50  1=1, MN 
IF<MaRKCI)>  49,5o,5o 

49  CF  =  F0CEA'N 
PP*F4t I ) 

IF<PP.LT. 101.35)  CF=FCDNT 
F4(I>  =  ,25*MY( I )*F5( I )*CF*F8( I ) 

Fl 0 ( I )=PP 

5 0  CONTINUE 


4  0 


STEP3P.71 
STEP3P.72 
STEP3P.73 
STEP3P ,74 
STEP3P , 75 
STEP3P.78 
STEP3P  . 77 
STEP3P, 78 
STEP3P,7? 
STEP3P.80 
STEP3P . 8i 
STEP3P.8? 
STEP3P.83 
STEP3P , 84 
STEP3P.6? 
STEP3P.86 
STEP3P.87 
STEp3p,88 
STEP3P.89 
STEP3p ,90 
STEP3P , 91 
STEP3P.92 
STEP3P ,93 
STEp3p , 94 
STEP3P , 98 
STEP3P ,9a 
STEP3P.97 
STEP3P.9S 
STEp3p,99 
STEP3P.1O0 
STEP3P , 1 0 1 
STEP3P , 10  2 
STEP3P , in  3 

•  STEP3P.104 
STEP3P , 1 05 
STEP3P ,106 
STEP3P.107 
STEP3P.108 
STEP3P.1P9 
STEP3P.110 
STEP3P , 1{1 
STEP3P.112 
STEP3P.113 
STEP3P.ll'4 
STEP3P.115 
STEP3P . 116 
STEP3p , 11 7 

STEP3P.H8 
APRl4,5l 
STEP3P , 120 

I 0N*»*»********»******STEP3P  lpl 

OR  OCEAN  SURFACE  STEP3P ,122 

Ps  IS  101.35  CB  OR  MOSTEP3P , ip3 
STEP3P ■ 124 
STEP3P,l25 
STEP3P.126 
STEP3P ,127 
S7EP3P . 1?8 
STEP3P , 129 
STEP3P , 130 
STEP3P. 131 
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STEP3P  (Continued) 


STEP3P ,132 
STEP3P.133 

.  .  STEP3P.134 

C* ***************************  SENSIBLE  HEAT***************** ****** *******STEP3P, 135 

C  HEaTINQ  PROM  OCEAN  SURFACE  WHEN  THE  air  IS  COLDER.  PQ  IS  the  TEMP  DIFFSTEP3P.136 
C  OCEAN  SURFACE  IS  ASSUMED  WHERE  STANDARD  PRESSURE  PS  IS  101.35  CB  OR  M0STEP3P.137 


CALL  RANWT(WS,F4) 
CALL  RANRD ( TS, F5 ) 


DO  59  1*2, Ml 

STEP3P .138 

Jl*  JM I N ( 1 ) * 1 

STEP3p,l39 

J21 JMAX(  I ) -1 

STEP3P.140 

K*( Jl-1)*M*I 

STEP3P . 141 

DO  59  Jb J1 i J2 

STEP3P ,142 

PP*F10 ( K ) 

STEP3p,l43 

PQ=F5<K)„H2*F2(K> 

STEP3P.144 

IF<PP.LT . 101 . 324 )  GO  TO  57 

STEP3P ■ 145 

IF<PQ,LE.0.°>  GO  TO  57 

STEP3P , 146 

PR  =  SQRTt ,25*MY(K)*( <F7(K*1)-F7<K-1> )*(F7(K*1)-F7(K-1) )♦ 

STEP3P , 147 

X  (F7(K*M)"F7(K"M> )*CF7(K*M)*F7(K»M) ) )  ) 

STEP3P.148 

F  5 ( K )  =  ,5E"2*h1*( ,1*PR-1. )*P0 

STEP3P ,149 

GO  TO  58 

STEP3P.150 

57 

F5(K)=0,0 

STEP3P ,151 

58 

K  =  K  +  M 

STEP3P ,152 

59 

CONTINUE 

STEP3P , 153 

r*******************#*******HUMIDITY  FQRECAST**************************STEP3P.154 

PM.1 

E1,E2  are  cOEFF  for  COMP  of  mean  STREAMFUNCTION 

STEP3P.155 

SSI 

,SS2,SS3  are  COEFF  FOR  COMP  OF  M£AN  DIVERGENCE.  D  IS  ZERO  OVER 

LANDSTEP3P.156 

THE 

HUMIDITY  IS  GIVEN  IN  ton  PER  SQUAREMFTER  AND  CENTIBAR 

STEP3P 1 157 

CALL  RANRJ(DIV1«F6> 

STEP3P . 158 

CALL  RANRD(DIV2#F7) 

STEP3P ,  159 

CALL  RANRD(HUM1#F8) 

STEP3P , 160 
STEP3P.161 

SSI  =  El+EM*C2 

STEP3P.1A2 

SS2  =  E2-EM*C1 

STEP3P , 163 

SS3  b.EM*C6 

STEP3P , 1*4 

DO  61  1=1, NN 

STEP3P.165 

PU  =  D 

STEP3P.166 

PP  =  FI  0 ( I > 

STFP3P.1A7 

IF<PP.LT. 101.35)  PO=0,0 

STEP3P.168 

F7(I)  =  F8(  I  )  * ( SS1*F6 (  T  )  *SS2*F7 (  I  )  +  F4  (  I  )  * ( PQ*SS3 )  ) 

STEP3P.169 

F  6  <  I  )  =>  E  M  *  F 1  (  I  )  ♦  E 1  *  F  2  (  I  )  ♦  F  2  *  F  3  (  I  ) 

STEP3P , I7O 

61 

continue 

STEP3P.171 

CALL  jACn9(F6#F0,FlQ,M,N> 

STEP3P.172 

DO  62  1=2* Ml 

STEP3P ■ 173 

J 1  =  J  M I N ( I  )*1 

STEP3P.174 

J2=JMAX( I ) -1 

STEP3P.175 

K  =<J1-1)*H*I 

STEP3P . 176 

DO  62  J  =  J 1 , J 2 

STEP3P.177 

F6(K)  =  F8<K+M>+F8<K-M)+r8<K*i>*F8<K-l>»4*F8(K> 

STEP3P.178 

<  =  K  +  M 

STEP3P.179 

62 

CONTIN  IE 

STEP3P. IPO 
5TEP3P.1S1 

PEPS  =  EPS*DELT 

STEP3P.182 

JJ  63  1  =  1,  'IN 

STEP3P , lfl3 

IF(MARK( I ) , GE . 0 )  GO  TO  63 

STEP3P , lfl4 

F8( I >»-DEPS*( #25*MY{ I)*F10( I >*F7( I)»ADlFT*My< I >*K6< I) ) 

STEP3P.135 

6  3 

CO  JTINS'E 

STEP3P.1«6 
STEP3P ,187 

CALL  RANRD ( HUM? i F 6 ) 

STEP3P.1P0 

00  65  1=1, MN 

STEP3p  ,189 

I F ( MARK ( I ) }  64,65,65 

STEP3P ,190 

64 

r 6  c I )  =  F 6  C  I  ) +F8 ( I ) 

STEP3P.191 

65 

CONTINUE 

STEP3P , I92 

C****»**********************»PRECIRITaT ION******«  ************  *********»*STPP3P. 193 


102 


STEP3P  (Continued) 


C  the 

C  PF.R 


66 


67 


6n 


69 

70 


71 


C  *  *  *  * 


1  79 
iSO 


3  37 


1  83 
18  ) 
190 


C 


rain  for  one  TIMESTEP  is  accumulated. 

SQUAREHETER 

DO  70  I  =  1 , M N 

I F ( M ARK ( I) )  66,70,70 

TEMP  =  IU*F2(I)*H4*F3(!) 

CALL  S AfilR  (  TEMP,  OSAT  ) 

PCI  =  F  6  (  I  )-,8*QSAT 
IF(PQ)  68,68,67 
F6(I>  =,rt*Q SAT 
F7(I}  =  .5F3*DELP*PQ 
IF<KT.EOil)  F7(I)=0,O 
I F  {  KT ,  EG .  2 )  F7(I)®2*F7(H 
OO  TO  70 
F7 ( I )  a.D 

PO  =  F6<  i >-,2*OSAT 
IF(PO)  69.70,70 
F 6 (  I  )  =  , 2*qs A  T 
-O-Nti^UE 

CALL  RAJ90<PREC# F8> 

JO  7 1  1  =  1,  IN 
F a <  I  )  =  F 8 ( I )  +  F7(  I  ) 
q  O  N  t  1 N  j  g 

CALL  RANJT(PREC,F8> 

CAuL  RAfHlHrf JM1#F8) 

CALL  RAinNT  <  H J M i ,  F 6  ) 

CALL  Ra  j>jT(HUM2,F8) 


THE  Rain  is  GIVEN  in  MM  OR  KSTEP3P,1q4 

STEP3P*Xa5 
STEP3P.196 
STER3P.197 
STEP3P  .198 
STER3P.199 
STFP3P.200 
STER3p,2nl 
STEP3P.202 
STEP3P,2q3 
STEP 3P  ,  2 0 4 
STEP3P.2D5 
STEP3p.2n6 
STEP3P.207 
STEP3P.208 
STEP3P , 2  0  9 
STEP3P.2l0 
STEP7P.211 
STEP3P . 21 2 
STEP3P . 21 3 
STEP3P,2i 4 

STEP3P  ,  21 5 
STEP3P.216 
STEP3P.217 
STEP3p.2i8 
STEP3P.219 


heat* 


,220 


CALL  RANROtDlVi, F6) 

CALL  RAN^j(DIV2,F8) 

JO  180  1=1, HN 
I F  (  -1ARK  {  I  )  )  179,180,18^ 

F  6  (  I  )  =  Tl*F4C  I  >*T2*F6U  )*T3*F8(  I  ) 
CONTIN  IE 

rJJ  190  1=1, MM 
l F ( M  A  R  K ( I ) )  I87,i9n,l90 


STEP3P.221 

STEP3P.2?2 

STEP 3P, 223 
3TEP3P.224 
STEP3P.225 
STEp3p. 226 
STEP3P.227 
STEP3P • 2?8 
STEP3P , 2?9 


RAIN  =  F  7 ( I ) /EPS/ JELT 


STEP3P.230 


IF  (  RA  I  I.lT.T'JL)  no  TO  188 
VERT  =  F6( I ) 

IF* VEHT.OT.-DELl)  GO  T^  18* 

)STAR  =  VERT 

IF*  VERf.GE. -CCS. AND.  VERT. LF, -DELI)  OSTAR  =  *-ABS(VERT*VERT/CC5) 
TEMP  a  H 6* F 2 ( j ) 

X  =  CC2**CC1-1./TEMP) 

E  =  E 0  *EXP ( X ) 

FSTAR  =  tE#TENP*E*(CC3-TFMp>/PMEAN/(PMFAN^TEHP*TtMP  *  CC4*E) 
HLAT  =  "  U*HL*OSTAR*FSTAR 
GO  TO  189 
HLAT  =0.0 

Fd*  I  )  =  F5< I >  *  HLAT 

continue 

CALL  RANRT(HEAT,F5) 

11=0 
12  =  0 

13  =  0 

14  =  0 


STEP3P.231 
STEP3P , 232 
STEP3P.233 
STEP3R.234 
STEP3P.235 
STEP^P.236 
STEp3p  t 237 

STEP3pt278 

STEP3P . 239 
STEP3P.240 
STEP3P.241 
STEP3P.242, 
STFP3P.243 
STEP3P.244 
STEP3P.245 
STEP3P.246 
STEP3P.247 
STEP3P.248 
STEP3P.249 
STEP3P . 230 


CALL  STEPEXT*F1,F2.F3,F4'jf5,F6,f7,F8,F9,F10,MY,MARK,M,N, 
X  n.12,13,14) 

C  ***** , 


A PR1 4 , 3b 
STFP3P.2S3 
>,  2S4 
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STEP3P  (Continued) 


CALL.  RANRD(J789,F10) 

UEPS  =  .25*DELT*EPS 
00  73  I s 1 j  MN 
I F ( MARK  < l ) )  72,73,73 

72  PQ  =  4*F< I ) / M Y ( I ) 

F3 ( I )  a-DEPS*(F10( I )  -  C**F  4 ( I ) * P Q ) 

F6(P  =  PQ*<A3#F4  <  1  >+Al*F5  <  P  > 

F7( I )=PQ*(B3*F4( I )+Bl*F5( I ) ) 

73  CONTINUE 

FALL  RaNRDC  J3,F10) 

CALL  RANRDC J12.F4) 

CALL  RANRU( J56#F5) 

JO  BO  1=1, MN 
I F  <  M ARK  C  I  )  )  79  »  8  0  #  80 
7o  PQ  =  F ( I  )  *F  < I  ) 

F6( I )  3-HEPS*  C  F  4 ( I )*PQ*(A2*F9( I).A1*F10( I ) )*F6( I) ) 
F7 (  I  )  a-l»EPS*(F5<  I  >*Pti*<Pjl*Fl  Q  f  I  >-B2*F9(  I  )  )»F7<  I  > ) 
80  CONTINUE 


STEP3P.285 
STEP3P.256 
STEP3P.257 
STFP3p,2s8 
STEp3p  t  259 
STEP3P.2A0 
STEP  3P . 26 1 
STEP3P , 26  2 
STEp3p,2f3 
STEP  3P . 264 
STFP3P . ?65 
STEP3P , 2*6 
STEP3P.267 
STEP3P.268 
STEP3P.269 
STEP3P.270 
STEP3P.271 
STPP3P.272 
STEP3P.273 
STEP3P.274 


0***»*************** ********* SOLUTION  OF  FORECAST  FQ, 
CALL  RANRD(HM3,p4) 

CALL  RANR0(Hl3,F5) 

CALL  RA Nrq ( H23  *  F 1 0 ) 

JEPSl=nELT*EPS 
DO  81  1=1* M  N 

F8< I )=F8( I >* HEPS 1/M Y( I  ) *F4 ( I ) 

F6( I )=F6< I )*UEPS1/MY( I >*F&<I) 

Ql  F 7 ( I ) =  F 7 ( I )*OEPSl/MY( I > *F*Q ( I ) 

CALL  RANRj(PSU2,F5) 

CALL  RANRD(PoI22#F10) 
n 

1)0  7  0  1  =  1,  IN 

89  F 2 ( I )  =  2*(F2( I )-F5( I ) ) 

Fj( I ) =  2  * ( F  3 ( I )-FlO< I ) > 

90  CONTINUE 


STEP3P.275 
STEP3P.276 
STEP3P,277 
STEP3P.278 
STEP3P.279 
STEP3P.2P0 
STEP3P.281 
STEP3P.2a2 
STEP3p , 2p3 
STEP3P.2P4 
STEP3P.2R5 
STEP3P ,  2  86 
STEP3P.287 
STEP  3p , 2gg 
STEP3P.2R9 
ST  EP3p , 2q  o 
STEP3P.291 


LAPEL  =  iOH  UELMSYS 
BTI1E  =  SECOND ( DU  4M Y ) 

NKMTE(6,8000)  RTIME 
BOOn  FORlATdH  .  *BTIME=  **  FlO.4) 

CALL  HEL  1SYS-(F2,F3,F6,  F7,My,F4#  Al,  A2.  B1,B2.  ALF ASYS.RFSSYS, 
X  ITSYS* M*N) 

DTI  IE  =  UTIME  -  SECOND t DUMMY ) 

*JRlTE(6,d005  )  LABEL*  DTI^E 
BOO'S  FORiATUH  ,  *TIME  to  EXECUTE  *«  AlO,  FlQ,4) 

IT(liKT)  =  I T  S  v  S 

JA‘LL  Ac>Ml)T(F2,F4,M,N,  ,S) 

CALL  ASMUT(F2,F4,M.N*-.5I 
CALL  ASM  NT  (  F  3 »  F  4  ,  M ,  N ,  ,*5) 
r ALL  ASMhT(F3,F4,M,U,-,5> 

c 

)0  100  I  =  1 1  M N 
I F (  1ARK( I ) )  99*100*100 

99  T F I L T  =  0 , 4 

I F <  1AR<( I ) .EQ.-10)  T F I L T  =  0 .7 
v  IF<MAR<< I ) .EQ.-l)  TFILTal, 

F 5(I)=F5(I)  *  TFILT*F2U) 

F10( I  )=FlO< I  )  ♦  TFILT*F3( I ) 

F4( I >  =  n/Myc I ) 

100  CONTINUE 
C 


APR  1 4 -57 
APR  14 ,58 
APR14 , 59 
APR1 4 , 6  0 
APR14 , 61 
STEP3P . 293 
APR14.62 
APRl4,ft3 
APRl^.64 
STEP3P.204 
JUN12.18 
JUN12.19 
JUN12 , 2  0 

JUN12,21 
STFP3P.295 
STEP3P.296 
STEP3P,297 
JUN12 .22 
JUN12.23 
JUN12 , 24 
JUN12.25 
JUNl? , 26 
STFP3Pt3nO 
STEP3P , 3nl 
STEP3P.302 
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STEP3P  (Continued) 


CALI  RANRDtPSIM2,F6) 

C 

DO  110 

109  Fl<l>  =  2*<F1<!’*F6(J))  «C2*F2*P  *C1*F3<P 

110  CONTINUE 

C 

CALL  HEL-1(Fl*F8«F4*  ALFAM’iRESM,  ITM,M,N) 

I T ( 2, KT)  «  I TM 
C 

CALL  ASMUT(Fl#F4fMJNf ,^) 

CALL  ASMUT(Fl-F4*M*N«*.5l 
00  120  1  a 1 j  N N 

IF<MARK< I ) )119.120,l2o 
119  TF I LTs  ,  4 

IFMARKI I ) .EG. -10)  TF  I  LTa  0 . 7 
I  F  <  “iARK  <  I  )  •  EQ  . -1  >  TFILT  =  1# 

F6(M*F6(I)  ♦  TFILT*FlU)  ♦  c2*f2(I) 

S-C  1*F3<  I  )) 

120  CONTINUE 

CALL  rANrD<STRM,F4) 

I F ( IVAR.EQ.O)  Q0  TO  156 
IF(KIND.NE.O)  GO  TD  156 
CALL  RANRD(ZMl,F7) 

CALL  RANRU(ZM2#F8) 

C 

RF  =  (KT-1)/ND 
FW  =  1,-WF 
3  =  9,806 
DO  130  1=1, MU 
I F <  MARK ( I ) )  129,  13  0, 1  3  n 
1^9  F 7 (  I  )  =  G*(FU*F7( I)+WF*F*( I >)/F( I ) 

130  C°NTINUE 

CALL  MIXF<F6#F7,MARK,WGTi , WGT2»WGT3#M,N) 
CALL  RANRD(STR1,F4) 

C 

CALL  RAnRQ(ZU,F7) 

CALL  RANRD(Z12'F6) 

00  140  I s 1 1 MN 
’ IF< lARK ( I ) )  1^9, l^o, l^n 
139  F  7 (  I  )  3  G*(Fri+F7(I)+wF*FA(I))/F(I) 

I4n  CONT  I  U  IE 

CALL  MIXF<F5#F7, MARK, WGTl# WGT2« WGT3.M,N> 
CALL  KANRD(STR2jF4) 

C 

CALL  RANRD(Z21*F7) 

CALL  RANRU<Z22*F8) 

no  I50  I s 1 , M N 

IFMARKU))  l49,l5o,l5n 
149  F7(I)  3  Fu*F7( I )^wF*Fa< I ) ) /F(  I ) 

150  CONTINUE 

CALL  i I XF { F 1 0 , F 7 , M ARK 1 Nq T1 , WG T2 , WQT 3 , M , N ) 

30  TO  156 

351  CALL  *1IXF(F6*F4.MARK.WGTi,WGT2iMGT3.M,N) 
CALL  RANRD<STRM,f4) 

CALL  I X  F  (  F  5  ,  F4 ,  M  ARK ,  WGTl ,  WGj2 ,  W3T3 ,  H N  ) 
CALL  RANRD(STR1*F4) 

CALL  MIXF(FlO,F4,MARK,WGTl.WGT2,WGT3,M.N) 
CALL  RANRD<STR2,F4) 

NER 


STEP3P.303 

STEP^P , 3  0  4 

STEP3P.305 

STEP3P , 3  n  6 

STEP3P.307 

STEP3P  ,  3 f! 8 

STEP3P , 3  0  9 

STEP3P.31 0 

STEP3P,3ll 

JU^12,27 

JUN12.28. 

STEP3P.312 

JUNl2,29 

JUn12,30 

JUN12.31 

JUNI12.32 

j  1  j  N 1 ? , 33 

JUNI^ , 34 

STEP3P.314 

FI  ELDS***  *************STEP3p.  31.5 

STEP3P,3l 6 

STEP3P.317 

STEP  3p , 31 8 

STEP3P.319 

STEP3P . 3?  0 

STEP3P,3^1 

STEp3p , 3?2 

STEP3P , 3?3 

STEP3P.324 

STEP3P , 3?5 

STEP3P.326 

JUN12  35 

STEP3P,3?8 

STEP3p.3?9 

STEP3P.330 

STEp3pt33i 
STEp3p.3l2 
STEP3P .333 
STEp3p,334 
STFp3pt335 
JLT12.36 
STEP3P.337 
STEP3P.338 
STEP3P . 339 
STEP3P , 3  4  0 
STEP3P.341 
STEP 3 P. 342 
STEp3p , 3  4  3 
STFp3pt344  , 
JlT12,37 
STEP3P.346 
STEP3P.347 
STEP3P ,348 
STEP3P.349 
STEP3P,350 
STEP3P.351 
STEP3P . 3^2 
STEP3P , 38  3 
STEP3p,3^4 
.355 
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STEP3P  (Continued) 


J5a  1 F  <  KT . LT . 3 )  00  TO  157 
CALL  FUNKDfPSIMl.Fl) 

CALL  RANPDiPSIll.F4) 

CALL  rt A N P d ( P S I  2 1 . F 7 ) 

157  CALL  RANWT  <  P5 IM1  j  F6 ) 

CALL  RA’UT(PSIUtF5) 

CALL  RA'JWT(PSl2l,F10) 

1F(KT.,LT.3>  00  TO  158 
CALL  RANRT ( PS  I M2<  Fl ) 

CALL  RA MT<PSJ12,F4> 

CALL  RAN. IT  <  PS  1 22,  F7  ) 

A-****************** 

call  hanr;j< j3,f!o) 
call  RA-JRD  ( HEAT,  F5) 

CALL  RA!WU<WS,F8) 

90  -  l./EPS/UELT 

DO  16Q  I = 1 • MM 

I r  <  MARK  <  I ) >  159,160,16' 

15?  TEFMl  =  F(I)*(PQ  *F2( I >♦ ,25*MYU ) *FlO < I > > -F5C I )*Si  *F8<I> 
T E ^ M 3  =  F(I)*(PQ  *F3 ( I  )  +  ,  2*5*MY< I  )  *F9  (I>>  -S2  *F8  < I > 

F2 ( I )  =  -  Al*TERMl  ♦  a2*TFRm2 

f6(I)  =  Hi* TERMl  -  B2*TESM? 


160  CONTINJE 


STEP3P.356 
STEP3P.357 
STEP3P . 358 
STEP3P,3^9 
STEP3P.360 
STEP3P.361 

STEP3P.362 
STEP3P ,3*3 
STEP3P.364 
STPP3P.3A5 
STEP3P.366 

}  OF  l)IVERGFNCE******************STFP3p(367 

STEP3P.368 
STEP3P.369 
STFP3P.370 
STEP3P , 37l 
STEP 3 P, 372 
STFP3P.373 
STFP3P.374 
STEP3P, 375 
STEP3P ,376 
STFP3P.377 
STFP3p.378 
STEP3P.379 
STFP3P.3A0 


CALL  R A WT(DIV1#F2> 
CALL  RANWT(D!V?#F3) 
CALL  R A  I R D ( P S I M 1 , F 1 ) 
CALL  RA-kJrd(PSI11(F5) 
CALL  RAMRJ(PSI?1*F10) 
C 

°r{HT  75 1 2 *  K T 
7512  FORMATUX,  3hKTs,  l3) 
170  cCiTlN'JE 


$TEP3P.3*1 

STEP3P.3R2 

STEP3p,3«3 

STEP3P,3«4 

STEP3P.385 

STEP3p.3fl6 

STEP3P.3R7 

STEP3R.388 

STEP3P.3R9 


RET  JRM 
END 


STEP3P.3Q0 
STEP3P , 391 
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SUBROUTINE  STEPEXT 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


a 


SUBROUTINE  STEPEXT  ( FI  *  F2*.  F3,  F4*  F5*  F6.  f7 #  F8,  F9*  FI  0  ,  MY  *  MARK #  H ,  N , 

1  11*12*13,14) 

THIS  SUBROUTINE  computes  the  CONTRIBUTION  from  the  HIGHER  ORDER  TERMS 
IN  THE  VORTICITY  EQUATION, 

ii  Indicates  the  vorticitt  advection  by  the  divergent  vind 

I  2  INDICATES  THE  RELATIVE  VORT I C I T Y*D I VERGENCE 

1 3  INDICATES  THE  VERTICAL  ADVECTION  OF  VORTIcITY 

14  INDICATES  the  twistingterm 
11=0  NO  CONTRIBUTION  u  DIFFERENT  from  0 


12=0  NO  CONTRIBUTION 
I  3  =  0  NO  CONTRIBUTION 
14=0  NO  CONTRIBUTION 


CONTRIBUTION  FROM  \% 

12  DIFFERENT  FROM  0  CONTRIBUTION  FROM  12 

13  DIFFERENT  FROM  0  CONTRIBUTION  FROM  13 
CONTRIBUTION  FROM  1 4 


14  DIFFERENT  FROM  0 

PSIM  IS  IN  FI#  PS  1 1  IS  IN  F?,PSI2  IS  IN  F3,WS  IS  IN  F4 
STEPEXT  NEEDS  10  FIELDS  IN  THE  FAST  CORE  MEMORY  F#MY  aNd  MaRK  MUST 
ALSO  BE  in  Fast  core  MEMORY 

DIMENSION  Fl<l>*F2(l) #F3fl) jF4<1)#F5(1)»F6(1)#F7<1)*F8(i) »T9(l) • 

1  F10<1).F(1),MARK(1),MY(1) 

common  f 

COMMON/COEFF/Al.  A2,A3,ni‘,B?,B3,Cl.C2,C3.C4,C5.C6,C7,C8,D.DFLP,EM. 

1  El*E2»Hl*H2*H3*H4*H5»H6*PMEAN*SlJS2'Tl*T2*T3»T4*T5* 

2  PO  #  PM  *  Pi 

COMMON/COEFF2/T6*T7,T8#T9*TlO*Tll,Tl2*Tl3#Tl4, 

X  Kl,K2,K34K4,K5<K6,K7,Kfl#K9,K10,Kll#Kl2,Kl3,Kl^*Kl5, 

X  K16,Kl7#KiA,Kl9,K20,K21#K22,K23,K24,K25,K26,K27,K28, 

X  K29*K30'm’31'K3?*K33»K34iK35<  K36#K37,K38 

COMMON/ECS/  PSIMliPSUl,oSl2l#PSIM2#PSI12*PSI22,HUMi,HUM2*DIVlJ 
20IV2<w:S,HEaTi J789* Jl2* J5f, j3,pS#TS,PREC#STRM,STRl,STp2,ZMl,Zll,Z2l 
3, ZM2, Z12,Z22,H13,H23,HM3JHM2, H12,H22,Hll,H2i; J4, VM, VI# V2 
REAL  MV 

REAL  Kl.K2iK3#K4#K5.K6#K7.K6#K9,KlO.KU#Kl2*Kl3#Kl4#Kl5. 

X  K16,K17,K1B#K19#K20#k2i,k2?,k23jk24#k25#k26,k27,k28# 

X  K29#K30#K3l#K32#K33,K34fK35#K36,K37,K3B 

<  I N  D  =  0 


APR14.67 
STEPEXT , 3 
STEPEXT, 4 
STEPEXT, 5 
STEPEXT , 4 
STEPEXT, 7 
STEPEXT, 8 
STEPEXf « 9 
STEPEXT, iO 
STEPEXT, H 
STEPEXT, 12 
STEPEXT, 13 
STEPEXT, 14 
STEPEXT, 15 
STEPEXT. 16 
STEPEXT . 1 7 
STEPEXT, 18 
APR14 ,68 
STEPEXT, 19 
STEPFXT , ?Q 
STEPEXT, 21 
STEPFXT, 72 
STEPEXT, *3 
STEPFXT, ?4 
STEPFXT, ?5 
APR14 ,69 
APR14,70 
APR14.71 
STEPFXT, 34 
STEPEXT, 35 
STEPFXT , 7  6 
STEPEXT, 37 
APR14 ,72 


MN=  1*N 

RESIDUE  =  ,  5E 4 
A  LF  A  =  1 , 4 

J  A  C  0  B  I  A  'I  J4  TO  SECONDARY  STO»  AGE  #  D  I  VERGED  I  k  S  TO  FAST  MEMORY 
CALL  RAMWT(J4,F9) 

CALL  RANHD ( D I VI #F5) 

CALL  RA-JRD ( D I V2 # F6 ) 
no  9  1=1* Mi 
IF( HaRK( I ) )  7,7*7 
7  F4U)  =  Q.O 
F5U  )  =  0  ,  0 
F6  < I ) =0 . 0 
9  CO  NT  IN  JE 

IF(KIN^.FU.O)  GO  TO  8 
CALL  B  I0VE(F4,M,N) 
call  R-IUVE(F5,M,N) 

CALL  BM0VE(F6,M, N) 

■8  CONjiNme 

I F( 14 ,EQ. 0  *  AND. 13 ,EQ, 0 . AND , I? ,E J. 0 , AND . 1 1 , EQ. 0 5  GO  TO  170 
I F  < I4.EQ.0)  GO  TO  44 
rO-IPJTE  THE  TW  I  ST  I NIGTERM#  I  4 
no  10  1=1, MN 

F7( I )sK3l*F5( I )*K3?*F6< I)*K33*F4<  I ) 
in  Fa<  I  )=K3<S*F5(  I  )+K37*F6(  I  )  +K38*F 4(1) 

CALL  GRADPR(F2#F7.F9#MARK#M,N) 

CALL  GRA!)PR<F3iF8#F10,MAPKjM,N> 

no  11  1=1, mN 

F9<  I )  =  0.5  +  MY( I ) *F9 (  I  ) 


STEPEXT.  38 
STEPFXT. 39 
STEPFXT, 40 
STEPFXT. 4I 
STEPFXT , 4  2 

STEPEXT, 43 

STEPEXT, 44 
STEPFXT. 45 
STEPFXT, 46 
STEPFXT. 47 
STEPEXT. 48 
STEPFXT. 49 
STEPEXT , 8  0 
STEPFXT. K1 
STEPEXT. S2 
STEPEXT. 53 

STEPFXT , 54 
STEPFXT. 55 
STFPFXT ,56 
STEPFXT , 57 
STEPFXT, 58 
STEPEXT  ,  *>9 
STEPFXT ,60 
ST FPFXT  #  ai 
STEPFXT , *2 
STEPFXT, 63 

STEPFXT, 64 
STEPFXT. 45 


s> 
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STEPEXT  (Continued) 


H  F 10< I  )  =  0  ,  5*MY ( I )*F10( I ) 

CALL  RANWT<H13.F9) 

CALL  RA.NWT  ( H23 * FlO  ) 
ju  20  1=1, MN 

F7< I )=«19*F5< I >+K20*F6( n+K21*F4C I > 

2 p  F8 ( I ) sK22*F5( I )+k23*f6( I)  +  K24*F4< I ) 

CALL  GRADPR(F2,F7,F9,Mark,m,N) 

CALL  GRA!JPR(F3-F8#FlO#MA9K*MiN) 

00  30  1=1, MN 

30  F7(  I  )=k25*f5(  I  >+k26*f6U)^k27*f4(  I  ) 

CALL  GRADPR(El,F7,FQ,M4R*<,M,N) 

)0  40  I=1*NN 

40  F8<I )=<Fd( I )+F9( I )+FlO( I ) )*Q ,5*MY( I ) 

CALL  RANUT(Hm3, F8) 

Li 0  TO  45 
44  46  1=1,  ^ 

46  F8(  I  )aQ  ,0 

CALL  RANWT(Hl3iF0) 

CALL  RA  iWT(H23-F8) 

CALL  RANRT(Mi3.F8) 

IFU^.FO.O.A'JD.  I3.FQ.O.A'JU.  Il.FJ.O)  GO  TO  I63 
SO  CALL  RFLV°R(Fi,F7,my,MARk,m, m> 

CALL  RFLVUR ( F2# F8* MY* MAR*' M* M ) 

CALL  RELV0R(F3*F9#MYjMARKiM,N) 

IFU2.EQ.O.A‘NL).  I3.FO.0)  GO  to  94 
I F ( I2.EQ.G)  GO  T°  65 

COMP JTATIO J  OF  THE  relative  vqrt ic i ty*d i vergence, 12 
JO  6  0  1=1,  MN 

Fl<  I  )=F5(  I  )*(Kl*F8(  I  )*'-<2*F9(  T  >+K3*F7<  I  >  )  +F  6  (  I  )  *  ( K4  * 
1  K6*F7< I ) )*F4< I >*<K7*F8<I  >+K6*F9( I )+K9*F7< I ) ) 

F2(I)=F5(I)*(K28*F8(I)*F7(I)>*F6(I)*K29*f6(I)+F4(I) 
6  0  F3( I >  =  F5(  I  )*K34*F9( I  )  +F6 (  I  > ★ (  K35*F  9< I )*F7< I) )♦ 

rOM FJTATIOI  IF  THE  VERTICAL  ADVECTION  OF  V0RTICITY#I3 
I  F  (  I  3  ,  E'J  ,  0  )  GO  TO  91 
JO  TO  HO 
65  .)U  70  1  =  1,  MN 
Fl (  I  >  =  0  .  O 
F  2  (  I )  e  0  .  »J 
70  F3(  I  >bo,0 
6n  )Q  90  1=1, MN 

Fl<  I  >sFj.(  I  )+F5<  I>*(KlO*Ffll  T  >*K11*F9<I  >*K12*F7<  I  >  )* 

1  F6(  I  )  * ( Kl3*F8 (  I  )*Kl4*F9(  I  )  *Rl5*F7 (  I  )  )  * 

2  F 4 ( I ) *  f  Kl6*F8 ( I >+Kl7*F9< I )*Kl8*F7< I ) ) 

F2(I)=F2(I>+F8(I)*<K31*FS(T)+K32*F6(I)*K33*F4(I)) 

9 H  F3( I )=F3( I J+F9U )#(K36*F5( I )^K37*F6( I )*K38*F4( I ) ) 

GO  TO  9I 
94  00  92  1=1* MN 
Fl(I)sO,0 
F2<  I  )  =  0,0 


STEPEXT, 66 
STEPFXT. 67 
STEPFXT ,68 
STEPEXT. 69 
STEPFXT. 70 
STEPE  XT , 7 1 
STEPFXT. 72 
STEPEXT, 73 
STEPFXT, 74 
STEPFXT , 75 
STEPEXT, 76 
STEPFXT. 77 
STEPEXT. 78 
STEPEXT. 79 
STEPFXT.«b 
STEPFXT.n 
STEPEXT, 82 
STEPEXT , «3 
STEPFXT, p4 
STEPFXT, 
STEPEXT. 86 
, A  P  R 1 4 , 73 
APR14 , 74 
APR14 ,75 
STEPFXT. 90 
STEPEXT ,91 
STEPFXT, Q2 
STEPEXT, 93 
F 8 (  I  )+K5*F9( I ) +STF  PF  XT ,94 
STEPEXT, q5 
*K3n*F8< I >  STEPEXT, 96 

F 4 (  !)*K30*F9< I)STEPFxt,97 
STEPFXT. g8 
STEPEXT. 99 
STEPEXT, 1 0O 
STEPEXT, 1 01 
STEPFXT, 1 0? 
STEPEXT, 1 03 
STEPFXT, 1 04 
STFPFXT. 1 05 
STEPFXT. 1 06 
STEPFXT. 1 07 
STEPEXT . 1 08 
STEPFXT. 109 
STEPF  XT , 1 10 
STEPFXT.1 11 
STFPFXT , 11? 
STEPFXT. 1 13 
STEPFXT. 1 14 


9?  F3 ( I ) =0 , 0 
9 1  CALL  RANHT(HM2,Fl) 

*  CALL  R A . H T ( H 1 2 , F 2 ) 

CALL  R A N ^ j ( H 2 2 j  P 3 ) 

i;  i,OMPiITATIUM  OF  THE  APVFCTInN  OF  VORTICITY  BY  THE  DIVERGENT  WIND  -It 

c  comp  j t g  forcingfhnction  fo~  thf  velocitypotential 

IF( Il.EO.O)  GO  TO  166 
95  JO  100  1=1, MN 

Fl C I }  =  <C^*F5(I )-cl*F6( I )^Cfi*F4(I ) ) / M  Y ( I ) 

F2( I )  =  F5( I ) / M  Y ( I ) 

ICO  F3<  I  )  =  F  6  < I >/MY< I ) 

C  SOLVE  THE  PO I SSONEOU AT  I  ON  0 Y  RELAXATION  IN  ORDER  TO  GET  VELOCITYPOT. 
CALL  RANRD<VM,F4) 


STEPFXT. 11,5 
STEPEXT. 116 
STEPFXT. 1 17 
STEPEXT, 118 
STEPFXT. 119 
STEPFXT , 1 20 
STEPF  XT  ,  t  2 1 
STEPFXT. 12? 
STEPFXT. 123 
STEPEXT. 124 
STEPEXT, 125 
STEPEXT, 126 
STFPFXT. 127 
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STEPEXT  (Continued) 


CALL  RANRD(Vl.FS) 

STEPEXT. 128 

CAUL  RANRUtV2,F6) 

STEPEXt . 129 

CALL  VELP0T(F4,Fl,M,N, RE 3 1  DUE,  ALFA ) 

STEPEXT, 130 

CAUL  VELP0T(F5.F2»M.N»RESIDUE,ALFA> 

STEPEXT. 131 

CALL  VELPOT < F6. F3, M.N, RES l DUE. ALFA > 

STEPEXT, 132 

CALL  RANWTCVM.F4) 

STEPEXt, 133 

caiIl  RA^kT ( V2, F 6 ) 

STEPEXT. 134 

CALL  RANWT ( VI » F5 ) 

STEPEXT, 135 

DO  HO  1=1,  MM 

STEPEXT. 136 

Fl<  I  )*F7C I )*F( ! >»2,*F8( I ) 

STEPEXT. l‘37 

110 

F2(  I  )=F7( I )*F( I )*2»F9( I ) 

STEPEXT. 1 38 

CALL  GRADPRtF4.F8.F3.M4RK, M.N) 

STEPEXT, 139 

CALL  GRA0PRtF5.Fl, FlO, MARK, M.N) 

STEPEXT. 140 

DO  120  1=1, MN 

STEPEXT, 141 

120 

Fl0<  I  )=-0 ,5*MY< I ) *<F3( I ) *FlO  < I ) ) 

STEPEXT. 1 42 

CALL  R4NWT ( Hll »FlO ) 

STEPEXT. 143 

CALL  GRADPRfF4.F9.F3. MARK, M.N) 

STEPEXT, 1 44 

CALL  GRA[)PR(F6.F2»F10,MaRK.M,N> 

STEPEXT. 1 45 

DO  130  1=1. MM 

STEPFXT.146 

130 

Fl0<  I)«-0,5*MY< I)*(F3(I)*F10(I>) 

STEPEXT. ijt7 

CALL  RANWT(H21,Fl0) 

STEPEXT. 14a 

DO  140  1=1, MN 

STEPEXT. 1 49 

Fit  I )=C3*{F7( I >+F( I ) ).C2*F8( 1 )+C4*F9f I >+C7*F  (  I  ) 

STtpEXT , 1 5  0 

F2(I>  =  -C2*<F7(I>*F<I>)+C’5*F8<I> 

STEPEXT.  1 5 1 

140 

F3( I )  =  C4* ( F7 ( I ) +F ( I ) )+C6*F9( I )+2.*C7*F<  I  ) 

STEPFXT.152 

CALL  GRAUPR(F4,f1,F7,MARK,M,N) 

STEPEXT. 153 

CALL  GRAUPRfF5.F2.F8, MARK, M.N) 

STEPFXT ■ 1 54 

CALL  GRAnPR(Fb,F3,F9,HARK,M,M) 

STEPEXT, 155 

DO  150  I =1  *  MV 

STEPEXT. 156 

n 

Fit  I  )  =  -0 ,5*HY(  I  )*(F7(  I  >+F8(  n*F.9(  I)  ) 

APR14 , 76 

CALL  HaNRD<H12.F2> 

STEPEXT. 15s 

CALL  RA 4RD<HM3,F3) 

STEPEXT , 1 5b 

CALL  RaNrd(H11»F4) 

STEPEXT. 160 

CALL  RANRD(H12*F5) 

STEPFXT. i6l 

CALL  RA'JRD(H13.F6) 

STEPEXT. 162 

CALL  RAmRU ( H21.F7 ) 

STEPEXT ■ 1 63 

CALL  RA,VRD(H22,F8) 

STEPEXT. 164 

CALL  RaNRD(H23.F9) 

STEPEXT, 165 

00  160  I=1»MN 

STEPEXT. 166 

Fit  I )  =  F1( I )*F2( I)*F3< J ) 

STEPEXT. 167 

F4( I )*F4t I >  *F5  C I )+F6< I ) 

STEPFXt . 1 68 

160 

F7 ( I )=F7( I )+F8( I >*F9< I  ) 

STEPEXT. 169 

GO  TO  190 

STEPFXT. 1 70 

3  63 

CALL  RANRDtHM3.Fl) 

STEPEXT. l7i 

CALL  RA NRD( Hl3» F 4 ) 

STEPFXT. 1 7? 

CALL  RANRD ( H23 . F7  ) 

STEPFXT. 173 

GO  TO  190 

STEPFXT. 1 74 

166 

CALL  RANRD ( HM2, F 2 ) 

STEPEXT. 175 

CALL  RANRU(HM3*F3) 

STEPEXT. 176 

CALL  HANRDt  H12, F5) 

STEPEXT. 177 

CALL  RAtJRD(Hl3,F6) 

STEPFXT. 1 78 

CALL  RAMRD(H22,F8) 

STEPEXT.179 

CALL  RANRD(H23*F9) 

STEPEXT.  1 80 

DO  167  1=1. MN 

STEPFXT. 181 

Fit  I )=F2( I )+F3( I ) 

STEPEXT, 1 82 

F 4 ( I ) =F5 ( I )+F6(  I  ) 

STFPEXT.183 

3  67 

F  7  t I ) =F8 ( I )+F9(  I  ) 

STEPFXT, 184 

•GO  TO  190 

STEPFXT, 1 85 

170 

DO  180  1=1. MN 

STEPEXT. 186 

Fit  I )  =  0  .  0 

STEPFXT. 187 

F4  t  I  )  =  0 , 0 

STEPEXT. 188 

1 8  0 

F7  (  I >  =  0  .  n 

STEPFXT. 189 
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STEPEXT  (Conti nued ) 


CALL  HA(\UT<HM3iFl> 

STEPEXT. 1 

CALL  RAN^T(Hl3,r4) 

STEPEXT. 1 

CALL  Ra’j,]T(H23,F7) 

STEPEXT. 1 

CALL  RANRU(PSIM1*F1) 

STEPEXT.  i 

CALL  RA.MRU(PSlU,F2) 

STEPEXT. 1 

CALL  RANrfD<PSl2l.F3) 

STEPEXT. 1 

CALL  RANr»D(WS,F4) 

STEPEXT. 1. 

CALL  RaNPU(HEaT,F5) 

STEPFXT.1 

CALL  RAi\lrt)J<  J4'F9> 

STEPEXT. i 

RETURN 

STEPEXT. 1 

£ND 

STEPEXT.  ? 

90 

91. 

92 

93 

94 

95 

96 

97 

98 

99 

'00 


